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ABSTRACT 

This  report  presents  the  relations  necessary  to  define  the  motion  of  a 
target  in  the  gravitational  field  of  the  earth.  In  order  to  express  this 
target  motion  in  the  frame  of  a radar,  ail  appropriate  set  of  coordinate 
systems  (and  transformations)  is  introduced.  Target  tracking  in  the  form  of 
a Maximum  Likelihood  Estimator  is  discussed.  The  problem  of  interceptor 
miss  distance  is  treated  from  the  standpoint  of  the  uncertainty  volumes 
associated  with  estimated  target  state  vectors. 
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SECTION  1 - - INTRODUCTION 


ARIES  is  a system  simulation  computer  program  developed  by  Lincoln 
Laboratory  to  study  radar  tracking  and  command  guided  intercepts  in  a 
realistic  radar  environment.  Written  in  FORTRAN  and  designed  for  execution 
on  the  CDC  6600  computer,  it  has  considerable  versatility  in  the  specifica- 
tion of  radar,  target,  tracking  and  environmental  models. 


ARIES  is  designed  to  be  a useful  analytical  tool  for  several  allied 
areas.  Originally  ARIES  was  used  in  the  strategic  RMD  area  to  estimate  the 
metric  state  vector  (position  and  velocity)  of  tracked  targets,  and  then  to 
extrapolate  ahead  in  time  to  determine  an  intercept  point  for  an  ICBM.  The 
radar  measurements  were  subject  to  environmental  effects  which  were  reflected 
in  the  intercept  miss  distances.  Refraction  and  scintillation  models  were 
used  in  ARIES,  and  the  effects  of  various  calibration  schemes  on  target 
location  accuracies  were  studied,  ARIES  was  also  used  to  study  the  problem 
of  multipath  in  low  angle  tracking  and  to  examine  the  effectiveness  of 
various  proposed  schemes  to  overcome  degradations  in  prediction  caused  by 
multipath.  The  use  of  ARIES  in  BMD  studies  was  terminated  in  the  summer  of 


More  recently,  a modified  endoatmospheric  version  of  ARIES,  known  as  the 
HWLPTR  (Hostile  Weapons  Location  Projectile  Tracking  Radar)  Program,  has  been 
used  in  the  tactical  area  for  the  evaluation  of  a radar’s  performance  in 
’’backtracking”  an  incoming  artillery  or  mortar  shell  to  determine  its  point 
of  launch.  The  resulting  point  estimation  error  CEP  values  also  assist  in 
the  evaluation  of  the  drag  model  error  effects  and  in  the  study  of  the  overall 
performance  of  hostile  weapon  location  systems. 

1.2  Program  Features 

Since  the  ARIES  Program  provides  a fairly  elaborate  system  simulation, 
it  is  useful  to  tabulate  the  various  features  incorporated  in  ARIES.  The 
major  components  of  the  simulation  are  summarized  below: 

1.  Target  trajectories ---accepts  input  of  target  state  vectors  in 
several  coordinate  systems.  Trajectories  with  launch  angle,  reentry 
angle  or  minimum  energy  constraints  from  a given  launch  location  to 
an  impact  point  are  also  available. 

2.  Radar  models — both  mechanically  steered  (dish)  and  phased  array 
radars  are  modeled.  Radar  sensitivity,  beamwidth,  frequency  and 
location  are  specified  by  inputs.  Range  and  angle  measurement 
precision  are  also  specified  by  inputs. 

3.  Radar  measurement  modeling 

a.  Target  modeling- --static  cross-section  measurements  on  real 

targets  are  used  in  conjunction  vrith  rigid  body  dynamics  (Euler’s 
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equations)  to  obtain  realistic  dynamic  RCS  simulations.  Constant 
and  sinusoidal,  as  well  as  an  analytic  cylinder,  RCS  options  are 
also  available. 

b.  Noise  and  propagation  effects — radar  measurements  are  corrupted 
by  receiver  noise  (S/N  dependent),  range- independent  noise  effects 
(jitter,  quantization,  etc.)  and  uncorrected  propagation  effects. 
Tropospheric  and  ionospheric  refraction  are  assumed  to  be  cor- 
rected to  within  random  percentages  (input  parameters)  of  the 
true  values.  Ionospheric  scintillation  and  multipath  effects 
corrupt  the  data  but  are  assumed  to  be  uncorrected  in  the 
measurement  model. 

4.  Trajectory  estimation  (Target  Tracking) — Maximum  Likelihood  Estima- 
tion (MLE)  of  the  target  trajectory  is  performed  based  on  measurement 
data  collected  at  specified  PRF’s  over  specified  track  intervals. 
Individual  measurements  are  weighted  according  to  their  measurement 
variances.  ARIES  could  be  easily  extended  to  use  recursive  tracking 
algorithms . 

5.  Target  Discrimination--- (Not  presently  implemented.)  Conceptually, 
discrimination  algorithms  would  be  implemented  to  determine  whether 
a particular  simulated  target  constituted  a threat  to  the  defended 
area. 

6.  Interceptor  modeling ---(Not  presently  implemented.)  Flight  charac- 
teristics of  one  or  more  interceptor  types  would  be  utilized  to 
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conduct  a command  guided  intercept.  Currently,  the  program  extra- 
polates the  estimated  and  true  target  state  vectors  to  various  time 
(or  altitude)  points  after  termination  of  track  to  obtain  miss 
distances.  Miss  distance  statistics  are  computed  from  the  accumu- 
lated miss  distances  observed  on  a series  of  Monte  Carlo  tests. 


In  addition  to  the  above  simulation  components,  ARIES  also  accommodates 
multiple  radars,  multiple  targets  and  multiple  track  intervals  on  a given 
simulation  run.  The  feature  of  making  many  Monte  Carlo  runs  for  a given 
scenario  permits  the  generation  of  meaningful  miss  distance  statistics.  A 
building  block/subroutine  program  structure  lends  itself  to  reasonable 
straightforward  modifications  of  or  additions  to  the  program. 

The  input/ output  of  the  ARIES  Program  is  engineer  oriented.  For  input, 
simulation  data  cards  are  conveniently  grouped  into  '‘packets"  (each  packet 
defines  a target  model,  a radar  model,  an  environmental  model,  etc.)  which 
the  engineer  may  simply  stack  up,  together  with  packets  specifying  the  desired 
simulation  "scenario".  For  output,  an  8V  x 11"  ARIES  Test  Report  (see 
Appendix  II  of  Reference  1)  is  generated  which  provides  the  engineer  with 
descriptions  of  his  input  model  parameters  and  scenarios,  along  with  the 
resultant  simulation  data  and  statistics.  The  outputs  are  ail  organized  into 
logical  sections  which  are  indexed  for  ready  reference.  Outputs  from  ARIES 
also  include  trajectory  plots  superimposed  on  a world  map,  plots  of  true  and 
measured  target  cross-section,  and  a radar  measurements  tape  containing  metric 
and  RCS  data  for  processing  by  other  programs. 
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1.3  Program  Documentation 

Th_  ARIES  Program  is  documented  in  three  separate  Lincoln  Laboratory 
Technical  Notes  as  follows: 

1.  The  ARIES  Program  - A General  Overview  and  Users'  Guide 

2.  The  ARIES  Program  - Coordinates,  Transformations,  Trajectories  and 
Tracking 

3.  The  ARIES  Program  - Analysis  and  Generation  of  Simulated  Radar 
Measurements 

The  first  report  presents  a general  discussion  of  the  ARIES  Program, 
including  the  logical  organization  of  the  program  and  descriptions  of  all 
subroutines.  All  of  the  options  available  to  a user  are  discussed  and  the 
methods  of  setting  up  the  input  "packets",  including  controls  to  activate  the 
various  options,  are  presented.  A typical  run,  including  a complete  ARIES 
Test  Report  (output),  is  discussed  in  this  first  volume. 

The  second  and  third  volumes  contain  all  of  the  relevant  mathematics  and 
the  models  used  in  ARIES.  Most  of  the  deterministic  mathematics  (coordinate 
systems  and  transformations,  trajectory  generation  and  estimation,  miss 
distance  calculations,  etc.)  are  in  the  second  report.  The  third  report  is 
primarily  concerned  withthe  generation  of  radar  measurements,  including  the 
corruptive  effects  of  noise,  radar  biases,  propagation  and  time-varying  radar 
cross-section. 


i] 
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1.4  Organization  or  This  Report 


In  Section  2,  some  aspects  of  the  WGS66  earth  model  are  presented.  This 
includes  a truncated  version  of  the  series  representation  of  the  gravitational 
field.  The  various  coordinate  systems  used  in  the  ARIES  Program  are  defined 
in  Section  3.  Target  accelerations  due  to  gravitational  forces  and  due  to 
atmospheric  drag  forces  (not  implemented  in  ARIES)  are  derived  in  Section  4. 

A simple  predictor -corrector  algorithm  for  the  integration  of  the  target 
equations  of  motion  is  also  derived.  In  Section  5,  all  of  the  transformations 
among  the  various  coordinate  systems  are  derived.  The  problem  of  defining 
realistic  target  trajectories  is  addressed  in  Section  6.  Only  targets  which 
are  launched  from  and  impact  on  the  earth  are  considered.  For  given  launch  and 
impact  points,  the  user  can  also  specify  the  launch  elevation  angle  or  the 
reentry  angle  or  he  can  request  a minimum  energy  trajectory.  The  trajectories 
are  first  derived  for  a central  force  field  (Keplerian  motion)  and  then 
perturbed  to  account  for  the  actual  gravitational  field.  Section  7 presents 
one  type  of  target  tracking  algorithm;  namely,  an  "after-the-fact"  maximum 
likelihood  estimate  of  the  trajectory  which  best  fits  the  accumulated  radar 
measurements.  From  the  state  vectors  generated  by  trajectory  fitting  on 
several  Monte  Carlo  runs  and  the  known  (true)  target  position,  one  can  develop 
an  error  covariance  matrix  representing  the  uncertainty  in  predicted  target 
position  at  an  assumed  intercept  point.  This  covariance  matrix  and  the  related 
handover  error  ellipsoid  and  error  sphere  are  presented  in  Section  8. 
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SECTION  2 --EARTH  MODEL 


The  model  for  the  earth  used  in  the  ARIES  Program  is  taken  directly  from 
the  1966  World  Geodetic  Survey  (WGS-66).  In  this  model  the  earth's  shape  is 
given  as  a surface  of  revolution  obtained  by  rotating  an  ellipse  around  its 
minor  axis.  The  resulting  surface  is  referred  to  as  an  ellipsoid  or  as  an 
oblate  spheroid.  As  part  of  this  same  sui/ey,  the  coefficients  required  for 
an  expansion  of  the  earth  gravitational  field  in  spherical  harmonics  were  also 
derived  from  the  measured  data. 

The  parameters  of  the  WGS-66  earth  model,  as  used  in  the  ARIES  Program, 
will  be  defined  and  summarized  in  the  following  sections. 


2.1  Earth  Ellipsoid 

The  representation  of  the  earth  as  an  ellipsoid  is  principally  a matter 
of  mathematical  convenience;  that  is,  actual  points  on  the  earth's  surface  will 
depart  in  varying  amounts  from  the  corresponding  points  on  the  ellipsoid. 
However,  on  average,  the  ellipsoid  will  be  a good  representation.  The  earth 
ellipsoid  is  defined  by  the  Equatorial  and  North  polar  radii  (semi -major  and 
semi -minor  axes).  These  values  are: 


Re  = 6378.145  Km 

(2-1) 

Rn  = 6356.760  Km 

(2-2) 

Two  other  parameters  of  the  ellipse,  the  flattening  factor  f and  the  eccentric- 
ity e,  are  also  of  interest 
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(2-3) 


*^1  - \^J  = yf (2-f)  = .081820 


(2-4) 


The  oblique  view  o£  the  earth  given  in  Figure  2.1  and  the  cross-sectional 
view  given  in  Figure  2.2  aid  in  the  definition  of  positions  on  the  earth 
surface.  Longitude  is  defined  as  the  angle  between  the  plane  containing  the 
Greenwich  meridian  and  a plane  containing  anv  other  meridian;  longitude  is 
measured  in  an  easterly  direction  from  the  Greenwich  meridian. 

The  earth  rotates  around  its  polar  axis,  as  indicated  in  Figure  2.1,  with 
an  angular  rate  of 


w = 7.29211515  x 10"  radians/ second 

V 


(2-5) 


There  are  two  latitude  angles  defined  for  the  earth  ellipsoid,  as  indicated 
in  Figure  2.2.  Both  angles  are  referenced  to  the  equatorial  plane.  The 
geodetic  latitude  <}>  is  defined  as  the  angle  between  the  local  normal  to  the 
earth  ellipsoid  and  the  equatorial  plane.  Mathematically,  it  is  defined  by 


1 

dz 

dn" 


{(l-f)y/ip,]j'1 


(2-6) 


(l‘f) ' 


ti/2  < < tt/ 2 


(2-7) 


Note  that  the  extension  o£  the  local  normal  to  the  equatorial  plane  does  not 
in  general  pass  through  the  center  of  the  earth  (exceptions  occur  when  <J>a0  and 
4»=tv/2)  . The  geodetic  latitude  is  the  usual  survey  latitude.  It  is  also  the 
latitude  determined  fran  astronomic  observations. 

Geocentric  latitude  <b  , on  the  other  hand,  is  simply  the  angle  formed  by 

v»* 

the  equatorial  plane  and  a vector  from  the  center  of  the  earth  to  a point  on 
the  earth  surface.  It  is  defined  by 


tan*  = - 
c n 


(2-8) 


or,  in  terms  of  geodetic  latitude, 


tan<j>c  = (1-f)2  tan<|> 


(2-9) 


The  radius  of  the  earth  at  an  arbitrary  location  on  the  surface  can  be 


expressed  as 


R = Jn^+z^ ' = 


ll-e2cos2 


(2-10) 


or,  in  terms  of  geodetic  la^i-uds, 


R = R_.  l-e2sin2<j>' 


(2-11) 


It  is  frequently  necessary  to  compute  the  altitude  of  a target  above 
the  earth  ellipsoid.  That  is,  for  a target  at  a location  (n0,zQ)  outside  the 


ellipse,  it  is  desired  to  find  the  minimum  of 


h - J(Vn)s  ♦ lVz)4' 


(2-12) 


subject  to  the  constraint  that  (n,z)  is  on  the  surface  of  the  ellipsoid: 


n2  + — = R2 


d-f): 


(2-13) 


A direct  attack  on  this  minimization  leads  to  a quartic  equation  which  is  very 
difficult  to  solve.  Fortunately,  the  eccentricity  of  the  earth  ellipsoid  is 
quite  small  and  a reasonably  accurate  approximation  can  be  made.  This  approx- 
imation is  simply  to  take  the  difference  between  the  distance  from  the  center 
of  the  earth  to  the  target  and  the  radius  of  the  earth  in  the  same  direction: 


h°  T°  yi-e2COS2A* 


(2-14) 


where 


rc  * fo  * zl 


(2-15) 


and 


COsA  = ~ 
*o 


(2-16) 


A detailed  evaluation  of  the  error  in  this  simple  relation  for  hQ  has  not  been 
made;  however,  sample  computations  indicated  worst -case  errors  of  the  order  of 
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10  meters  for  altitudes  around  2000  JOn  and  about  4 meters  for  altitudes  around 
700  Km. 

For  the  purposes  of  the  ARIES  Program,  the  above  model  for  hQ  is  probably 
adequate.  However,  a refinement  of  the  relation  can  be  made  by  taking  first 
order  perturbations  of  n and  z around  hQ.  Such  a procedure  yields 


hi  = 


(l-e2COS2A) 

^l-(2-e2)e2cos2Ak 


(2-17) 


Sample  calculations  indicate  wcrst-case  errors  of  only  3 meters  for  h = 2000  Km 
and  only  0.5  meters  for  h = 700  Km. 

Assorted  coordinate  systems  associated  with  the  earth  model  and  appropriate 
transformations  between  each  pair  of  these  coordinate  systems  will  be  defined 
in  Sections  3 and  5. 


2.2  Earth  Gravitational  Field 

The  gravitational  potential  of  the  earth  is  generally  expanded  in  a 
series  of  spherical  harmonics.  The  most  important  component  of  the  potential 
is  the  point  mass  potential 

VQ(r)  = | (2-18) 

where  the  gravitational  constant  G is 

G = 3.986012  x 1014  (meters) 3/ (second) 2 (2-19) 
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and  r is  the  radius  from  the  center  of  the  earth.  Additional  contributions  to 
the  gravitational  potential  result  from  the  mass  distribution  of  the  earth. 
These  perturbations  of  the  point  mass  potential  are  expanded  in  spherical 
harmonics.  One  series  of  harmonics  depends  only  on  latitude;  these  are  called 
the  zonal  harmonics.  All  remaining  terms  of  the  harmonic  series  have  a 
dependence  on  both  latitude  and  longitude;  these  are  called  tesseral  harmonics. 

The  ARIES  Program  only  considers  the  zonal  harmonics  in  the  gravitational 
potential  model.  That  is,  the  potential  is  given  by 


(2-20) 


til 

where  Cn  is  the  coefficient  or  the  n — zonal  harmonic  and  Pn(0  is  the 
Legendre  polynomial  of  the  first  kind  of  order  n.  The  angle  A is  the  angle 
between  the  equatorial  plane  and  the  vector  from  the  center  of  the  earth  to 
che  target  position.  The  first  order  harmonic  is  zero  as  a result  of  the 
symmetry  of  the  gravitational  field.  Since  the  coefficient  of  the  second 
zonal  harmonic  is  roughly  three  orders  of  magnitude  greater  than  all  other 
coefficients,  the  ARIES  Program  is  further  simplified  to  include  only  the 
point  mass  potential  and  the  second  zoual  gravitational  harmonic.  (Additional 
zonal  harmonic  coefficients,  through  the  ninth,  are  available  In  the  ARIES 
Program;  they  are  easily  implemented  in  Subroutine  GRAVTY.) 
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Values  for  the  harmonic  coefficients  in  the  gravitational  potential  model 


C2  * 1.08271  x 1C"' 
C3  = -2.630  x 10'6 
= -2.35  x 10'6 
C5  = -0.265  x 10“6 
Cg  = 0.66  x 10"6 

C7  = -0.46  x 10  6 
Cb  = +0.53  x 10'6 
Cg  = 0.24^  x 10'6 


(2-21) 
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SECTION  3- -COORDINATE  SYSTEMS 

Several  different  coordinate  systems  are  used  in  the  ARIES  Program.  The 
various  operations  performed  iu  the  Program  (e.g.,  trajectory  integration, 
target  tumbling,  radar  measurement  generation)  all  have  preferred  coordinate 
frames  for  their  evaluation.  In  the  sections  below  the  various  coordinate 
systems  employed  by  the  ARIES  Program  are  defined.  The  transformations 
between  the  various  coordinate  systems  are  presented  in  Section  5. 

3.1  Earth  Centered  Inertial  Frame  (ECI) 


A primary  inertial  frame  employed  by  the  ARIES  Program  is  the  Earth 
Centered  Inertial  (ECI)  frame.  The  origin  of  this  frame  is  at  the  center  of 
the  earth,  the  z-axis  passes  through  the  North  Pole,  and  the  equatorial  plane 
of  the  earth  lies  in  the  x-y  plane.  At  the  reference  time  (typically  the 
launch  time  of  the  target) , the  x-axis  of  the  ECI  frame  passes  through  the 
Greenwich  meridian.  The  x,  y,  and  z axes  shown  in  Fig.  2.1  form  an  ECI  frame 
if  the  axes  are  held  fixed  in  space  rather  than  rotating  with  the  earth. 

The  equations  of  motion  for  ballistic  targets  have  their  simplest  form  in 
the  ECI  frame.  Consequently,  it  is  the  preferred  frame  for  trajectory  genera- 
tion and  trajectory  integration. 
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3.2  Earth  Centered  Fixed  Frame  (ECF) 

This  frame  is  essentially  the  same  as  the  ECI  frame,  except  that  it 
rotates  with  the  earth.  Usually  the  x-axis  is  defined  to  point  through  the 
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equator  at  the  Greenwich  meridian  and  the  y-axis  is  defined  to  point  through 
the  equator  at  90°  East  longitude. 

The  ECF  frame  is  a convenient  frame  for  the  definition  of  locations  on 
the  earth  surface.  Radar  sites,  target  launch  points  and  target  impact  points 
are  fixed  points  in  an  ECF  frame. 

3.3  Radar  Coordinate  Systems 

There  are  four  different  coordinate  systems  associated  with  the  radars. 

Two  of  these  are  peculiar  to  a phased  array  radar  (XRp,  Y^p,  Zp^  and  R,  U,  V) . 
One  coordinate  system  (R,  A,  E)  is  associated  with  a conventional,  mechanically 
steered  radar.  The  fourth  coordinate  system  is  an  earth  surface  fixed  frame 
or  radar  XYZ  frame.  Radar  measurements  are  typically  made  in  either  RAE  or 
RUV  coordinates.  The  cartesian  frames  are  useful  for  many  computations  and 
provide  the  necessary  cartesian  frames  for  definition  of  RAE  and  RUV  coordin- 
ates . 


3.3.1  Radar  Cartesian  Coordinates  (XYZ) 

The  Radar  Cartesian  Frame  is  defined  such  that  the  X-Y  plane  is  tangent 
to  the  earth  ellipsoid  and  the  Z-axis  is  pointed  outward  along  the  local 
vertical.  The  X-axis  points  East  and  the  Y-axis  points  North.  This  set  of 
coordinates  is  illustrated  in  Fig.  3.1.  Usually  the  origin  of  this  coordinate 
frame  will  coincide  with  the  radar  location  or  with  the  target  launch  point. 

If  the  radar  is  not  located  on  the  earth  ellipsoid,  but  is  actually  at  height 
H above  the  ellipsoid,  then  the  radar  X-Y  piano-  will  be  parallrl  to  the  tangent 
plane  to  the  earth  ellipsoid. 
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3.3.2  Radar  RAE  Coordinates 


The  RAE  coordinates  are  defined  relative  to  the  radar  cartesian  coordin- 
ates as  indicated  in  Fig.  3.2.  R is  the  distance  (range)  from  the  origin  (the 
radar)  to  the  target.  The  azimuth,  A,  is  the  angle  from  the  Y-axis  (North)  to 
the  projection  of  the  range  vector,  R,  into  the  X-Y  plane.  Azimuth  is  measured 
clockwise  from  North  toward  East.  The  elevation  angle,  E,  is  measured  (positive 
up)  from  the  X-Y  plane  to  the  range  vector,  R. 

3.3.3  Phased  Array  Radar  Face  Cartesian  Coordinates 

The  cartesian  coordinates  (X^p,  Y^,  Z^p)  are  defined  with  respect  to  the 
phased  array  radar  face  as  indicated  in  Fig.  3.3.  The  radar  face  lies  in  the 
Xrp'Yrf  plane,  with  the  XRp-axis  horizontal  and  the  Y^-axis  pointing  upward. 


The  Zpp-axis  points  outward  along  the  normal  to  the  array  face. 

The  orientation  of  the  phased  array  face  is  defined  by  the  azimuth  and 
the  elevation  of  the  phased  array  boresight  (i.e.,  the  phased  array  Z-axis). 
Azimuth  and  elevation,  in  this  case,  are  measured  in  the  same  sense  as  defined 
above  for  radar  RAE  coordinates. 

3.3.4  Phased  Array  RUV  and  Rag  Coordinates 

The  RUV  and  Rag  coordinate  systems  are  def  i relative  to  the  Phased 
Array  Radar  Face  Cartesian  Frame.  R is  the  dii,’  .ice  (range)  from  the  origin 


(the  radar)  to  the  target.  Target  direction  in  angle  is  measured  by  the 
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cosines  are  designated  U and  V,  respectively.  Note  that  conical  surfaces 
about  the  XRp  and  about  the  Y^p  axes  are  defined  by  U = constant  and  by  V = 
constant;  the  intersection  of  these  cones  defines  the  direction  to  the  target. 
(In  the  usual  spherical  coordinate  system  (R04>)  0 defines  a cone  about  the  z- 
axis  and  <f>  defines  a plane  containing  the  z-axis;  the  intersection  of  the 
conical  surface  and  the  plane  defines  target  direction.) 

The  angles  a and  3 can  also  be  used  to  measure  the  direction  to  a target 
with  respect  to  the  Phased  Array  Radar  Face  Cartesian  Frame.  In  this  case,  a 
and  6 are  the  arcsines  of  U and  V,  respectively.  Note  that  (tt/2-o)  and  (tt/2-3) 
are  the  direction  cosine  angles  with  respect  to  the  XRF  and  YRp  axes.  The 
phased  array  broadside  direction  is  characterized  by  U = V = 0,  which  corres- 
ponds to  direction  cosine  angles  of  ir/2  (and,  therefore,  0 = 3*0). 


3.4  Target  Coordinate  Systems 

The  target  tumbling  and  RCS  generation  programs  utilize  several  different 
coordinate  systems  which,  are  unique  to  this  particular  subject.  In  particular, 
the  following  frames  are  employed: 


1.  Principal  body  axes  frame  - this  frame  has  its  origin  at  the  center 
of  mass  of  the  object  and  rotates  with  the  object. 

2.  Momentum  frame  - this  is  an  inertial  frame  with  its  z-axis  aligned 
with  the  momentum  vector.  (The  momentum  vector  is  invariant  in  an 
inertial  frame.) 
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3.  Local  horizontal  frame  - this  frame  has  its  z axis  in  the  direction 
of  the  ECI  velocity  vector  and  its  y-axis  normal  to  the  trajectory 
plane. 

These  frames  are  inter-related  by  the  target  Euler  angles.  The  local 
horizontal  frame  is  directly  defined  from  the  target  position  and  the  target 
velocity  in  the  ECI  frame. 

Complete  details  about  these  coordinate  systems  and  transformations  among 
them  are  presented  in  Section  7 of  Reference  2. 


SECTION  4 --BALLISTIC  EQUATIONS  OF  NOTION 

In  this  Section  the  equations  o£  motion  for  a target  on  a ballistic 
trajectory  are  presented.  The  first  part  gives  the  required  relationship  for 
target  motion  in  the  presence  of  only  the  gravitational  field  of  the  earth. 

Atmospheric  drag  effects  are  not  included  in  the  acceleration  model  for 
the  ARIES  Program,  but  it  seems  desirable  to  explain  how  they  might  be  added 
to  the  model.  This  is  done  in  Section  4.2. 

The  last  part  (4.3)  explains  the  predictor-corrector  technique  for 
trajectory  integration  in  the  ARIES  Program. 


4.1  Target  Acceleration  Due  to  Gravity 

A target  located  at  a position  (x,y,z)  outside  the  surface  of  the  earth 
experiences  gravitational  forces  which  result  in  the  following  target  accelera- 
tions* 


x = W(r) 


where  V is  the  gradient  operator 


V ~ Sc  9x  + *y  3y  + *z  dz 


(4-1) 


(4-2) 


and  V(r)  is  the  gravitational  potential  function  defined  in  Section  2.2.  Tar- 
get position  is  given  in  the  ECI  coordinate  frame.  The  three  cartesian 
components  of  the  acceleration  vector  can  be  shown  to  be: 


*Lower  case  script  characters  denote  vector  (or  column  matrix)  quantities. 
Upper  case  script  characters  denote  matrices. 


X = 

Gx  „ 
• — 

(4-3) 

r3 

y = 

-^81 

(4-4) 

r3 

z = 

Gz  „ G „ 

' *T  gl  ' “7  s? 

rJ 

(4-5) 

where 

81  “ * - S (r)  CnPn*l(slnA) 

n=z  \ / 

9 /r  \n 

82  = £ (/)  (^(sinA)  (4-7) 

n=2  \ / 


In  these  relations  P^C  •)  is  the  derivative  of  the  Legendre  polynomial  with 
respect  to  its  argument.  Also, 

r = y/x2+y2+z2’  (4-8) 

and 

sirA  = (4-9) 

The  following  recursion  relation  for  the  derivatives  of  the  Legendre  poly- 
nomials 

PiUl(a)  = (2  + H}[aPA(a)  ' Pn-l(a)]  + PA-l(n°  (4'10) 
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is  useful  for  the  computation  of  and  g2.  Note  also  that 
Pq(cO  - 0,  P'(ci)  = 1 and  PJ(a)  = 3a. 


(4-11) 


Given  the  target  position  and  the  target  velocity  in  the  ECI  frame  at  a 
particular  reference  time,  together  with  the  above  definitions  for  target 
accelerations,  the  motion  of  the  target  in  the  ECI  frame  is  completely  defined 
over  the  entire  trajectory.  This  statement,  of  course,  ignores  any  effects 
due  to  atmospheric  drag. 

Subroutine  GRAVTY  is  used  to  compute  the  gravitational  acceleration 
components  of  a target  state  vector. 


4.2  Target  Acceleration  Due  to  Atmospheric  Drag 

The  ARIES  Program  does  not  include  the  target  deceleration  effects  caused 
by  atmospheric  drag.  However,  for  completeness,  a brief  discussion  of 
atmospheric  drag  is  included. 

Atmospheric  drag  results  in  a component  of  acceleration  in  the  direction 
of  the  velocity  vector  given  by: 


*D  = ~ 1 Cj)  A p v2/(w/g) 


(4-12) 


where 


Cp  = coefficient  of  drag  of  the  reentry  body 
A = cross-sectional  area  of  the  reentry  body 
p = atmospheric  density 


r 


v = body  velocity  (magnitude) 
w = body  weight 
g = acceleration  of  gravity 

The  parameters  Cp,  A and  w depend  on  the  shape  and  size  of  the  target  and 
would  have  to  be  inputs  to  the  ARIES  Program.  Atmospheric  density  is  a 
function  of  several  parameters  including  altitude,  temperature , latitude,  and 
longitude.  An  appropriate  atmospheric,  density  profile  would  also  have  to  be 
input  to  the  program,  either  as  an  analytic  formula  or  as  a table.  Body 
velocity  is  in  either  an  ECF  frame  or  an  earth  surface  fixed  (i.e.,  radar  XYZ) 
frame.  With  the  following  definition 


%•  Vv 


target  acceleration  due  to  drag  can  be  written  as 


(4-13) 


h>x 

aDy 

aDz 


-CDV* 
"CD  vy 
-Ci>vz 


(4-14) 


where  v , v and  v are  the  three  cartesian  components  of  velocity  in  an  ECF 

a y 

frame.  These  accelerations  must  then  be  transformed  to  the  ECI  frame  and 
added  to  those  due  to  gravity  to  obtain  total  target  acceleration  in  the 
atmosphere.  (The  atmosphere  is  typically  defined  to  cover  the  altitude  regime 
from  0 to  300  Kft  (91.44  Km).  Re-entry  altitude  is  defined  to  be  300  Kft  in 
the  ARIES  Program.) 
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The  addition  of  atmospheric  drag  would  result  in  a better  modeling  of 
target  motion.  However,  it  would  add  considerable  complexity  to  the  ARIES 
Program,  particularly  in  the  trajectory  generation  and  trajectory  estimation 
areas,  without  any  significant  impact  on  the  radar  measurement  modeling  which 
is  a major  part  of  the  ARIES  Program. 

4.3  Predictor -Corrector  Trajectory  Integration 

Trajectory  integration  in  the  ARIES  Program  is  performed  with  a predictor- 
corrector  integration  scheme.  The  procedure  for  integration  from  tx  = 0 to 
t2  - At  is  explained  below. 

If  the  accelerations  at  the  beginning,  a(0) , and  end,  a(At),  of  the  inter- 
val were  known,  the  acceleration  over  the  interval  could  be  approximated  by  a 
straight  line: 

a(t)  = a(0)  + -a-^'a(0-)  t (4-15) 

The  position  and  velocity  at  time  At  are  then  obtained  from 

At  t 

x(At)  = x(0)  + v(Q) At  + / dt  / a(x)dx  (4-16) 

0 0 

= x(0)  + v(0)  At  + a(0)(At)z/2 
+ [a(At)-a(0)](At)2/6  (4-17) 

and 


v(At)  = v(0)  + a(0)At  + [a(At) -a(Q) ]At/2 


(4-18) 
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In  practice,  of  course,  a (At)  is  not  known  so  the  integration  must  be 
done  in  two  steps.  First,  based  on  the  known  target  position  at  t = tx,  the 
accelerations  ax,  a and  a2  are  computed  from  the  gravity  model.  Then,  target 
position  and  velocity  are  predicted  for  time  t " t2  - tj  + At: 


Xp(tx+At)  = x(tx)  + v(tj)At  + a(t1)(At)2/2 
VpCti+At)  = vCtj)  + a(tj)At 


(4-19) 


(4-20) 


From  the  predicted  position  at  t * t2,  the  end-point  accelerations  are  computed. 
Then  a final  correction  of  the  predicted  (integrated)  position  and  velocity 
is  made: 


x(t2)  = + [a(t2)-a(t1)](At)?/6 

v(t2)  = V (t2)  + [a(t2) -3(^)1  At/2 


(4-21) 

(4-22) 


This  procedure  is,  of  course,  performed  simultaneously  on  all  three  cartesian 
coordinates  of  the  ECI  state  vector. 

The  extrapolation  (integration)  of  a target  state  vector  from  tl  to  t:,  is 
performed  by  Subroutine  EXTRAP  (X,T2).  The  components  cf  the  input  target 
state  vector  X are: 


X(l)  = Tl,  valid  time  of  the  state  vector 
X(2-4)  = ECI  positions  in  x,  y,  z 

X(5-7)  = ECI  velocities  in  x,  y,  z 
X(8-10)  = ECI  accelerations  in  x,  y,  z 
X ( 11)  = target  altitude 
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Integration  is  performed  with  a maximum  step  size  of  1 second  to  prevent  a 
buildup  of  errors.  (Tests  have  indicated  that  At  = 1 second  gives  very 
accurate  results.)  The  target  state  vector  is  valid  for  the  time  T2  when  it 
is  returned  to  the  calling  program  from  Subroutine  EXT FIAT. 


SECTION  5- -COORDINATE  TRANSFORMATIONS 


The  ARIES  Program  utilizes  several  different  coordinate  systems  as  defined 
in  Section  3.  Results  of  operations  performed  in  one  coordinate  system  are 
frequently  required  in  a different  coordinate  system.  For  example,  a trajec- 
tory is  integrated  in  an  ECI  frame,  but  it  is  used  in  a Radar  RAF  (or  RUV) 
frame  for  the  generation  of  radar  measurements. 

This  Section  presents  all  coordinate  transformations  relevant  to  trajec- 
tory generation,  trajectory  estimation  and  tracking  error  evaluation.  Some 
additional  transformations  involving  body  coordinate  systems  for  RCS  generation 
are  given  in  Ref.  2,  Section  7. 

5.1  Between  ECI  and  ECF 

At  the  reference  time  (usually  t = 0) , the  ECI  and  ECp  frames  are 
coincident.  Both  have  their  x-axes  through  the  equator  at  the  Greenwich 
meridian  and  their  z-axes  through  the  North  Pole.  At  a later  time  t the  ECI 
frame  is  unchanged  but  the  ECF  frame  has  rotated  with  the  earth  through  an 
angle  u>  t.  Let  the  subscript  I indicate  a component  in  the  ECI  frame  and  the 

v 

subscript  F a component  in  the  ECF  frame.  Now  it  is  easily  shown  that  the 
rotation  around  the  z-axis  by  an  amount  w t results  in  the  transformation 


XF 
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or,  by  defining  the  coefficient  matrix  as  G(t)  and  the  position  vectors  as  x^ 

and  xr:* 

F 

xF  = GCt)Xj  (5-2) 


If  one  takes  the  time  derivative  of  Eqn.  (5-1),  the  following  transformation 
of  the  velocity  vector  is  obtained: 


sF  = sctjjtj 


* 

0 

0) 

e 

0 

’Wi” 

"we 

0 

0 

*1 

= G(t) 

^rwexi 

0 

0 

0 

«■ 

h 

(5-3) 


where  the  dot  over  a vector  (or,  a vector  component)  indicates  a derivative 
with  respect  to  time.  The  second  term  in  Eqn.  (5-3)  results  from  the  contin- 
uous rotation  of  the  ECF  frame  relative  to  the  inertial  frame.  For  the  ARIES 
Program  we  require  only  position  and  velocity  transformations  between  ECI  and 
ECF  frames.  Target  accelerations  are  used  only  for  trajectory  integration 
(ECI  frame)  and  therefore  do  not  have  to  be  transformed. 

The  transformation  matrix  G(t)  represents  an  orthonormal  transformation. 
Consequently,  G(t)  is  orthogonal;  i.e.,  the  inverse  is  equal  to  the  transpose: 


G(t)  = G(t)T 


The  reverse  transformation,  from  ECF  to  ECI,  is  therefore  given  by 


*j  = G(t)Txp  (5-4) 


*Script  characters  denote  matrices  (lower  case  for  column  matrices;  upper  case 
for  N x M matrices) . 
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for  positions  and  by 


*1  - 0(t)T  j*F 


w 0 

e 


-*e  0 0 xp  =G(t)‘  >VVp  (5‘5) 


‘F'Vf 


for  velocities. 

Subroutine  ECIECF  (XI,XF,TIME)  is  used  to  transform  the  state  vector  XI 
in  EC I coordinates  to  the  state  vector  XF  in  ECF  coordinates.  TIME  is  the 
interval  in  seconds  since  the  two  frames  were  coincident.  Subroutine  ECFECI(XE, 
XI, TIME)  performs  the  inverse  transformation. 


5.2  Between  ECF  and  Radar  XYZ 

Consider  a radar  located  at  geodetic  latitude  <|>,  longitude  X (measured 
east  from  the  Greenwich  meridian)  and  at  height  h above  the  earth  ellipsoid. 
The  process  of  transforming  from  the  ECF  frame  to  a radar  XYZ  (East,  North, 

Up)  frame  can  be  visualized  as  the  following  sequence  of  operations: 

1.  Rotation  by  angle  A about  the  z-axis. 

2.  Rotation  by  angle  $ about  the  y-axis. 

3.  Translation  by  the  local  earth  radius  along  the  x-axis. 

4.  Rotation  by  angle  (4> ~4>c)  about  the  y-axis. 

5.  Translation  by  h along  the  x-axis. 

6.  Coordinate  interchange  to  obtain  the  desired  East,  North,  Up  sequence, 


The  appropriate  matrices  for  these  operations  are: 


1. 

cosx 

sinX 

8 = -sinX 

cosX 

0 

0 

2. 

cos<j>c 

0 

C = 0 

1 

;sin*c 

0 

(5-6) 


(5-7) 


= I 0 > where  Rg 


(l-e2cos2 


cos(4>-4>_)  0 sin(<j>-<j>  ) 


- sin(4>  _4>c)  0 cos(<j>-<|>c) 
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If  the  subscript  R refers  to  the  radar  frame,  then  we  can  write 


xR  = {EPC8}Xp  + E{k+V*e} 


(5-12) 


When  the  indicated  matrix  multiplications  are  performed,  the  resulting 
relations  are: 


XR  ” + k 


(5-13) 


where 


A = 


-sinA  cos A 0 

-sin<i>  cosA  -sin<$>  sinA  cos<j> 

cos<j)  cos  A cos<}>  sinA  sin<j> 


(5-14) 


and 


b = 


Rgh  sin(<|>-(J>c) 
-Reh  cos(*-«c)-h 


(5-15) 


Since  b simply  represents  a translation,  it  does  not  enter  into  the 
velocity  transformation.  For  the  velocity  transformation  f.'om  ECF  to  radar 
XYZ,  we  have  simply: 


*r  - AXp 


(5-16) 


The  matrix  A is  orthogonal,  so  the  inverse  transformation,  Radar  XYZ  to 
ECF,  is 


and 


.T  jT, 
Xp  = A Xp-A  b 


(5-17) 


*F  ” ^ *R 


(5-18) 


The  transformation  matrices  A and  b are  initially  computed  in  Subroutine 
SET(JP(W,A)  where  W is  a three  vector  containing  radar  longitude,  geodetic 
latitude  and  height  above  the  earth  ellipsoid.  The  A and  b matrices  are 
returned  in  the  transfer  vector  A.  This  vector  A is  then  passed  to  Subroutines 
ECFXYZ(XF,XR,A)  and  XYZECF(XR,XF,A)  to  perform  the  actual  coordinate  trans- 
formations . 

5.3  Between  Radar  XYZ  and  RAE 

The  Radar  Range,  Azimuth  and  Elevation  are  defined  in  Section  3.3.2.  They 
are  related  to  the  Radar  XYZ  coordinates  as  follows: 

R = yjx2+y2+z2 

A = tan  ^x/y),  0 ^ A < 2n  (5-19) 

E = sin  1 (z/R),  -V2  - E < 


The  proper  quadrant  must  be  determined  for  A so  that  A ranges  from  0 to  2tt, 
Time  derivatives  of  R,  A and  E yield: 
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(5-20) 


These  transformations  from  Radar  XYZ  to  Radar  RAE  are  performed  by  Subroutine 
XYZRAE. 

The  transformation  from  RAE  to  XYZ  is  performed  by  Subroutine  RAE XYZ  with 
the  use  of  the  following  relations: 


x = R sinA  cosE 


y = R cosA  cosE 
z = R sinE 


(5-21) 


x = S sinA  + yA 
y = S cosA  - xA 
i = R sinE  + RE  cosE 


(5-22) 


where 


S = ft  cosE  - RE  sinE 


(5-25) 


5.4  Between  Radar  XYZ  and  Phased  Array  XYZ 

The  definition  for  the  Phased  Array  XYZ  coordinate  system  was  given 
previously  in  Section  3.3.3.  Basically,  the  orientation  of  the  array  face  is 
defined  by  the  azimuth  and  the  elevation  1;  of  the  normal  to  the  array  face. 
The  array  face  lies  in  the  x.p-yp  plane  with  xp  horizontal  and  yp  up.  The  zp 
axis  is  aligned  with  the  array  normal.  (The  subscript  p has  the  same  meaning 
as  the  subscript  RF  used  in  Sections  3.3.3  and  3.3.4.) 

If  we  start  with  the  array  face  coordinates  coincident  with  the  radar  XYZ 
frame  (that  is,  the  array  is  pointed  up),  then  the  following  sequence  of 
transformations  gives  the  desired  Radar  XYZ  to  Phased  Arrav  XYZ  transformation: 


-1  0 0 1 0 


0 0 lJLo  cosF.p  sinEv 
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0-10  0 sinEp  -cosEp  sinAp  cosAp  0 yR  (5-24) 
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sinAp  sinEp 
sinAp  cosF.p 


-cosA  sinE.  cosF. . 


cos  A,.  cosE„  sinE, 
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Note  that  when  Ap  = Ep  = 0 (i.e. , the  phased  array  is  pointed  north),  Xp  = -xR, 
Yp  = zR  and  Zp  - yR,  as  would  be  expected  from  the  definition  of  the  phased 
array  XYZ  coordinate  system. 

The  velocity  transformation  is  simply 


*p  = Arp*R 


(5-27) 


Since  the  r trix  is  orthogonal,  the  transformation  matrix  for  con- 
verting from  Phased  Array  XYZ  to  Radar  XYZ  is  simply  the  transposed  matrix  A^1 
Subroutine  XYZXRF(XR,XP,Y)  is  used  to  convert  a state  vector  xR  in  Radar 
XYZ  to  a state  vector  Xp  in  Phased  Array  XYZ  coordinates.  The  reverse  trans- 
formation is  performed  by  Subroutine  XRFXYZ(XP,XR,Y) . In  both  Subroutines  the 
two  element  vector  Y contains  the  azimuth  Ap  and  elevation  Ep  of  the  phased 
array  normal. 


5.5  Between  Phased  Array  XYZ  and  RUV 


The  Phased  Array  range  and  "angle5'  coordinates  RUV  are  defined  by: 


R = Ji c2+y2+z2' 
V P 7P  p 


U = Xp/R  = sina,  -1  < U s 1 


V = y /R  = sinB,  -1  s V s 1 


(5-28) 


The  velocity  relations  are  obtained  by  taking  derivatives  of  Eqn.  (5-28): 


* BMStEiaihiV;;; ft  uirLi-  v.- t at  j,  •. ^ 


ft  = VpVpVp)/r 


u = OyUR)/R 


V = ap-VR)/R 


(5-29) 


The  transformation  from  Phased  Array  XYZ  coordinates  to  Phased  Array  RUV 
coordinates  is  performed  by  Subroutine  XRFRUV(XP,RUV) . 

The  inverse  transformation,  from  RUV  to  Phased  Array  XYZ,  has  the 
following  relations: 

X = RU 
Y = RV 

Z = r/i-U2-V2  = RW 


and 


(5-30) 


X = RU  + RU 
Y = RV  + RV 

2 - fav  - «(™±m 


(5-31) 


This  transformation  is  performed  by  Subroutine  RUVXRF(RUV,XP) 


5.6  Target  Coordinate  Systems 

The  target  tumbling  and  RCS  generation  programs  utilize  several  different 
coordinate  systems  which  are  unique  to  this  particular  subject.  These  coordi- 
nate systems  are  briefly  defined  in  Section  3.4.  A complete  description  of 
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these  coordinate  systems  and  the  transformations  among  them  are  presented  in 
Section  7 of  Reference  2. 

5.7  Other  Combinations 

The  coordinate  transformations  defined  above  consider  only  transfonnations 
between  quite  closely  related  coordinate  systems.  Other  transformations  may 
be  obtained  by  the  use  of  a sequence  of  the  above  transformations.  For  example, 
to  transform  a target  state  vector  in  the  EC I frame  to  a target  state  vector 
in  the  Phased  Array  RUV  frame,  the  following  sequence  of  transformations  would 
be  performed : 

1.  EC I to  ECF  (ECIECF) 

2.  ECF  to  Radar  XYZ  (ECFXYZ) 

3.  Radar  XYZ  to  Phased  Array  XYZ  (XYZXRF) 

‘ Phased  '"ray  XYZ  to  RUV  (XRFRUV) 

The  nr.mes  in  parentheses  are  the  names  of  the  ARIES  Subroutines  which  perform 
the  indicated  transformations.  This  particular  transformation  (ECI  to  Phased 
Array  RUV)  and  the  transformation  from  ECI  to  Radar  RAE,  are  qui.e  important 
and,  consequently,  they  are  available  by  calling  a single  Subroutine  ECIRAD 
(XI  ,W,N,TYPE) . I1''  is  t!  input  state  vector  in  ECI  coordinates,  W is  the  out- 

put state  vector  in  Radar  coordinates,  N is  the  identification  number  of  the 
Radar  (for  use  in  multiple  Radar  coverage)  and  TYPE  identifies  the  radar  as 
either  a dish  radar  or  a ;ed  array  radar. 
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The  coordinate  transformations  defined  in  Sections  5.1  to  5.5  provide  all 
the  information  required  to  transform  a state  vector  from  one  frame  to  any 
other  frame.  More  than  or. a subroutine  may  have  to  be  called,  as  in  the  above 
example,  but  the  key  point  is  that  the  capability  does  exist. 


5.8  Radar  Coverage  Limits 

A conventional  dish  radar  usually  has  the  capability  of  illuminating  and 
tracking  a target  over  the  entire  hemisphere  above  the  local  horizon  plane. 

On  the  other  hand,  a phased  array  is  usually  a permanently  fixed  installation 
(the  array  normal  direction  does  not  change).  Consequently,  the  coverage  of  a 
phased  array  radar  face  is  limited  to  a fraction  of  the  hemisphere.  Typical 
phased  array  coverage  lies  in  the  range  of  1/4  to  1/3  of  a hemisphere. 

For  a phased  array  radar,  in  particular,  it  is  important  to  test  the 
target  state  vector  to  determine  whether  the  target  lies  within  the  assumed 
field  of  view  of  the  phased  array.  This  test  is  performed  on  the  RAE  target 
state  vector  in  Subroutine  LIMITS (N, RAE, FLAG) . The  radar  number  N and  the  tar- 
get state  vector  in  RAE  coordinates  are  inputs  to  the  subroutine.  FLAG=1  is 
returned  if  the  target  is  in  the  field  of  view;  FLAG=0,  otherwise.  For  a dish 
radar,  the  test  simply  amounts  to  checking  that  target  elevation  is  greater 
then  zero.  A phased  array  radar,  on  the  other  hand,  is  assumed  to  have  the 
capability  of  viewing  a sector  in  azimuth  and  elevation  about  the  bores ight 
pointing  direction.  For  example,  if  the  boresight  azimuth  of  the  phased  array 
is  Ap,  then  one- third  hemisphere  coverage  is  achieved  by  applying  azimuth 
bounds  of  Ap  - tt/3  and  Ap  + tt/3  and  elevation  bounds  of  0 and  it/ 2.  Both 


azimuth  and  elevation  are  checked  against  the  phased  arrav  radar  coverage 
bounds  to  determine  whether  the  target  is  within  the  field  of  view.  Note  that 
it  is  easier  to  define  phased  array  coverage  limits  in  azimuth  and  elevation, 
rather  than  the  normal  phased  array  U and  V coordinates. 
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SECTION  6- -TRAJECTORY  GENERATION 

The  ARIES  Program  requires  the  capability  of  generating  a wide  range  of 
target  trajectories  to  simulate  various  radar/target  geometries . For  the 
present  requirements,  it  is  sufficient  to  consider  only  target  trajectories 
launched  from  and  impacting  on  the  surface  of  the  earth.  In  addition, 
atmospheric  drag  is  ignored  and  the  trajectories  take  the  short  way  around  the 
earth  from  launch  to  impact. 

The  trajectory  generation  process  requires  the  following  inputs: 

1.  Launch  longitude  and  geodetic  latitude 

2.  Impact  longitude  and  geodetic  latitude 

3.  Type  of  trajectory 

a.  Launch  angle  specified,  or 

b.  Re-entry  angle  specified,  or 

c . Minimum  energy 

The  procedure  used  in  the  generation  of  a trajectory  is  to  first  get  a good 
trajectory  estimate  using  Kepler's  equations  for  satellite  motion  in  a central 
force  field.  Perturbations  of  this  first  estimate  are  then  made  to  account 
for  the  zonal  harmonics  of  the  gravitational  field  and  to  account  for  the 
oblateness  of  the  earth  ellipsoid. 

Section  6.1  contains  a brief  summary  of  the  equations  which  define  the 
motion  of  an  object  in  an  elliptic  orbit  around  the  earth  (central  force  field 
only).  This  is  followed  by  a discussion  of  initialization  of  the  trajectory 
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state  vector  subject  to  the  launch,  impact  and  type  of  trajectory  constraints 
in  Section  6.2.  Finally,  the  procedures  employed  in  perturbing  the  initial 
state  vector  to  obtain  the  desired  precise  state  vector  are  presented  in 
Section  6.3. 

6 . 1 Equations  for  an  Elliptic  Orbit 

In  this  Section,  a set  of  equations  is  developed  for  target  motion  in  an 
elliptic  orbit  which  intercepts  the  surface  of  the  earth.  Figure  6.1  illustrates 
the  geometry  of  the  problem  as  viewed  in  the  plane  of  the  trajectory.  (Note 
that  this  plane  is  only  defined  in  the  ECI  frame.)  Since  two  of  the  trajectory 
types  considered  involve  specification  of  angles  (launch  angle  or  re-entry 
angle) , the  following  derivations  of  the  orbit  parameters  contain  launch 
elevation  as  an  explicit  variable. 

Most  standard  physics  textbooks  contain  the  derivation  of  the  orbital 
motion  of  an  object  in  the  presence  of  the  point  mass  gravitational  potential. 
They  are  summarized  here  for  the  gravitational  potential  of  the  earth. 

Radius  (center  of  earth  to  ellipse): 


= ajl-e2). 
l-ecos0 


(6-1) 


where  a is  the  semi-major  axis,  e is  the  eccentricity  of  the  orbit  and  o is  the 
central  angle  measured  from  apogee  (see  Fig.  6.1).  Note  that  the  true  anomaly 
is  measured  from  perigee  and  is  equal  to  (n-0).  The  radius  can  also  be  defined 
in  terms  of  the  eccentric  anomaly  Ec: 
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r = a(l-ecosEc) 


(6-2) 


If  these  two  definitions  for  r are  equated,  it  can  be  shown  that  true  anomaly 
and  eccentric  anomaly  are  related  by: 


tan(Ec/2)  = Jirf'tan(V^) 


(6-5) 


Two  other  parameters  are  also  of  interest: 


Mean  motion: 


n - \Jc,/a3' 


(6-4) 


Mean  anomaly: 


M * n(t-t) 


(6-5) 


M * E -esinE„ 
c c 


(6-6) 


where  t is  the  time  of  passing  perigee  and  G is  the  gravitational  constant  de- 
fined by  Eqn.  (2-19).  For  a full  orbit,  Ec  ranges  from  0 to  2n.  M also  goes 
from  0 to  2it,  so  the  orbital  period  is 


T = 2-n/n 


(6-7) 


The  total  energy  and  the  angular  momentum  are  constants  for  the  orbit.  That 


^m(f-2+r?-Q2)  - = constant 


(6-8) 
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and 


r20  = H = constant 


(6-9) 


Now  we  want  to  express  the  various  orbital  parameters  in  terms  of  the 
launch  elevation.  A first  step  is  to  demonstrate  that  the  angle  a in  Figure 
6.1  is  equal  to  2E.  As  a means  of  proving  this  statement,  write  the  equations 
for  the  circle  and  ellipse  as 


x2+y2  = R2 
1 e 

(6-10) 

(x-ae)2  + — = a2 

1 -o2 

(6-11) 

Then  compute  the  unit  vectors  tangent  to  the  circle  and  ellipse  at  the  launch 
point  (Recos0Q,Resin0o) : 


f,  - ^xsinoo  - XyCos0o;  circle 


(6-12) 


'ScSin0O  ‘ ^y(C0SVe) 
^e 


yi+e2-2ecos©  1 


; ellipse 


(6-13) 


The  elevation  angle  E is  now  determined  as  the  inverse  cosine  of  the  dot 
product  of  the  unit  vectors  i and  -t  : 


1-ecosO, 


cosE  = = 

c e 


i/l+e2-2ecosO  1 
v o 


(6-14) 


46 


Use  of  the  identity  cos2E  = 2cos2E-l,  yields  | 

i-e2-2ecos0  (l-ecos0  ) 1 

cos2E  = 2 (6-15)  | 

l+e2-2ecos0^  1 

0 1 

Application  of  the  Law  of  Cosines  to  the  triangle  formed  by  the  vectors  from 
the  foci  of  the  ellipse  to  the  launch  point  (see  Figure  6.1)  yields: 

4a2e2  = R2  + (2a-RJ2-2R  (2a-RJcoso  (6-16) 

6 C G 6 

This  relation  can  be  rewritten  in  the  form 

l-e2-2ecos0„(l-ecos0  ) 

cosa  = ^ 21  (6-17) 

l+e2-2ecos0„ 
o 

when  the  definition 

a = R„(l-eco.s0  )/(l-e2)  (6-18) 

V U 

is  used.  (This  is  Eqn.  6-1  evaluated  at  the  launch  point  r = Re,  0 = 0 .) 

Inspection  of  the  expressions  for  cosa  and  cos2E  demonstrates  that  a = 2E. 

The  Law  of  Sines  can  also  be  applied  to  the  triangle  to  obtain  an 
alternate  expression  for  the  major-axis  of  the  ellipse: 


2a-RQ  R 

e _ e 

sinOQ  ~ sin(TT-2E-G0) 


(6-19) 
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(6-20) 


Also  from  this  triangle  we  have  that 


2ae  2a"Re 


(6-21) 


The  use  o£  Eqn.  (6-20)  in  Eqn.  (6-21)  yields,  after  some  manipulation,  the 


following  relation  for  eccentricity: 


(6-22) 


If  the  total  energy  at  apogee  is  equated  to  the  total  energy  at  perigee,  then 
it  can  be  shown  that  the  angular  momentum,  H,  is  given  by 


H = \Z^a(l-e2)' 


(6-23) 


Now,  the  total  energy  at  launch  can  be  equated  to  the  total  energy  at  apogee 
(r  = a(l+e)  at  apogee) : 


1 2 _ G_  = 1 H2 

Re  2 a2(l+e)2 


(6-24) 


G /2aRe 

PT  V — 5 


(6-25) 


The  use  of  Eqn.  (6-20)  for  a gives  the  desired  result: 


- B r sii 

VL  ]/  Re  [sin©0+s: 


(6-26) 


The  eccentric  anomaly  at  launch  (see  Eqn.  6-2)  is  given  by 


ecl  " cos 


-iK" 

at 


and  the  eccentric  anomaly  at  impact  is 


(6-27) 


ECI  = 2tt_ECL 


(6-28) 


From  and  Eq,  the  mean  anomalies  at  launch  and  impact  can  be  computed. 


The  difference  between  the  two  is 


AM  = Mj-M^  = 2(Tr-Er^+esinE(-.jj) 


(6-29) 


Use  of  this  value  of  AM,  then  yields  an  expression  for  time  of  flight 


At  = tj-tL  = AM/r.  =v/^“ 


(6-30) 


This  completes  the  summary  of  the  elliptic  orbit  parameters  for  object 
motion  in  a central  force  field. 


6.2  Trajectory  Initialization 

The  summary  presented  above  in  Section  6.1  gives  the  parameters  of  a 
Keplerian  elliptic  orbit  in  a convenient  form  for  use  in  obtaining  an  initial 
estimate  cf  the  target  state  vector.  In  this  Section,  the  equations  will  be 
specialized  to  accommodate  particular  trajectory  constraints. 

Let  the  launch  point  have  longitude  X^  and  geodetic  latitude  4^  and  let 
the  impact  point  have  values  Xj  and  <j> j . The  geocentric  latitude  is  readily 
computed  from  the  geodetic  latitude.  A unit  vector  from  the  center  of  the 
earth  to  the  launch  point  is  given  by 

- <x  cos<t>cL  cosX^+ty  cos<|>cL  sinXj+-t2  sin4>cL  (6-31) 

During  the  flight  of  the  target  the  earth  rotates  an  amount  u t^,  where  uj  is 

e r e 

the  rotation  rate  and  t £ is  the  time  of  flight  (initially  unknown) . The  impact 
longitude  is  rotated  by  this  amount 

Xj , - X^  + (°gt£  (6-32) 

and  a unit  vector  in  the  direction  of  this  rotated  impact  point  is 

t = lx  cos<|>cI  cosXj.^y  cos4>cI  sinXj,+^z  sin4>cI  (6-3 ’X 

The  earth  central  angle,  2oq,  between  the  launch  and  impact  points  is  given 
by 


cos2e„  = i, 
o I L 


(6-34) 
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or 

cos2go  = cos4>cj  cos^  cos(X^-Xj  ,)+sin4>cj  sin4>cL  (6-35) 

It  is  j'.-  eful  to  have  a vector  in  the  trajectory  plane  which  is  orthogonal 
to  the  launch  unit  vector: 

n = (-t.  x tp)  x ^ = lv  - cos20Q^  (6-36) 

This  vector  is  tangent  to  the  unit  sphere  and  in  a direction  from  launch 
toward  impact.  The  magnitude  of  this  vector  is  readily  shown  to  be  |sin(20Q) | , 
so 

in  - yj-q-  3 Uj'-cos20o^L)/sin2Go  (6-37) 

Since  the  trajectories  are  assumed  to  always  be  the  short  way  around  the  earth, 
0 < eQ  < tt/2  and  the  magnitude  sign  on  sin2©Q  is  unnecessary. 

The  velocity  vector  at  launch  lies  in  the  plane  defined  by  the  unit 
vectors  and  -cn>  If  the  launch  speed  is  v^  and  the  launch  elevation  is  E, 
then  the  launch  velocity  vector  in  the  ECI  frame  can  be  written  as 

VL  = viJ(*Lsin^+'SicosE)  (6-38) 

center  of  the 
at  the  launch 
By  means  of  the 
vector  can  be  put 


This  velocity  vector  combined  with  a vector,  R^-c^,  from  the 
earth  to  the  launch  point  gives  a complete  ECI  state  vector 
point  (Rl  is  the  radius  of  the  earth  at  the  launch  point) . 
coordinate  transformations  defined  in  Section  5,  this  state 
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into  any  desired  coordinate  system.  The  remainir j part  of  the  problem  -- 
namely,  assigning  appropriate  values  to  and  E--is  considered  below. 

6.2.1  Minimum  Energy  Option 

It  is  frequently  desirable  to  fly  the  target  along  the  trajectory  which 
requires  the  minimum  amount  of  energy  at  launch.  This  is  achieved  by 
minimizing  the  launch  speed.  The  mininum  launch  speed  in  the  Eel  frame  occurs 
for  a launch  elevation  of 

E - J - ^ (6-39) 

Unfortunately,  this  is  minimum  energy  in  the  ECI  frame,  rather  than  an  earth 
surface  fixed  frame  where  the  minimum  likely  will  occur  under  different 
conditions  due  to  coriolis  effects.  If  the  components  of  the  launch  position 
are  denoted  by  x^>  y^,  and  the  components  of  the  vector  by  nx,  n^,  nz, 
then  the  velocity  components  in  the  ECF  frame  are 

vx  = vL[sinE(xL/RL)  + cosE  nx]  + u>eyL 
Vy  = vL[sinE(yL/RL)  + cosE  ny)  - u>exL 

vz  = VjJsinE(Zj/RL)  + cosE  n2l  (6-40) 

The  magnitude  squared  of  the  velocity  in  the  ECF  frame  is 

v2  “ \ + <4(x l+y20  + 2u)evL  cosE^Lnx‘xLny)  t6'41) 
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This  relation  indicates  that  the  elevation  angle  which  minimizes  v2  will  only 
coincide  with  the  elevation  which  minimizes  v^  when  the  last  term  is  zero. 

The  last  term  contains  the  z-component  of  the  cross  product  of  and  i which 
is  non-zero  unless  the  trajectory  plane  (ECI)  contains  the  North  Pole. 

The  procedure  employed  for  minimization  of  v2  is  to  search  for  the  eleva- 
tion angle  E which  results  in  zero  slope: 


U.  V ~ 

W = '2vl 


Jvlcos(2E+0q)  + 

| sxnOQ  + sin(2 


BcosEcosg, 


(6-42) 


where 


6 ■ “e(yLnX'\ny) 


Note  that  the  minimum  occurs  at  E 5 tr/4  - 0Q/2  for  a Polar  orbit  (g=0).  An 
initial  guess  of  E * tt/4  - 0 /2  is  made,  followed  by  steps  of  0.01  radians 
until  the  zero  of  dv2/dE  is  bracketed.  Then  a two-point  linear  approximation 
of  dv2/dE  is  used  in  an  iterative  scheme  to  locate  the  value  of  E which 
minimizes  v2.  This  value  of  E is  then  used  to  compute  the  orbit  parameters  a, 
e,  E^,  AM  and  At.  The  new  time  of  flight  At  is  compared  to  the  previous 
value  and  the  whole  procedure  repeated,  unless  the  difference  is  less  than  one 
second. 

Finally,  the  ECI  state  vector  with  the  position  vector,  R^,  and  the 
velocity  vector  vL  [given  by  Eqn.  (6-38)]  is  transformed  to  an  earth  surface 
fixed  XYZ  frame  (Radar  XYZ).  The  final  optimization  of  the  state  vector  (which, 
among  other  things,  accounts  for  flattening  of  the  earth  and  for  the  effects 
of  the  second  gravitational  harmonic)  is  discussed  in  Section  6.3. 
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6.2.2  Launch  Angle  Specified 

This  option  requires  the  determination  of  an  initial  velocity  vector  vL 
in  the  ECI  frame  which  when  transformed  to  an  earth  surface  fixed  frame  (at 
the  launch  point)  has  the  proper  elevation  of  the  velocity  vector.  That  is, 
if  vv,  v and  v,  are  the  components  of  the  velocity  vector  in  the  earth 

x y u 

surface  frame,  then 


Ef  = 


(6-43) 


must  be  made  equal  to  the  desired  launch  angle,  E^. 

An  iterative  procedure  for  the  determination  of  the  velocity  vector 
follows.  For  a given  value  of  the  elevation  E in  the  ECI  frame,  compute  the 
corresponding  launch  speed  (Eqn.  6-26),  and  launch  velocity  vector  (Eqn.  6-38). 
Convert  the  ECI  launch  state  vector  thus  determined  to  an  earth  surface  XYZ 
frame  at  the  launch  point.  Then  compute  the  launch  elevation  angle  from 
Eqn.  (6-43)  and  compare  it  to  the  desired  value  E^.  If  Ep  is  not  close  enough 
to  Ed  (e.g.,  0.001  radians),  choose  a new  value  of  E and  repeat  the  calcula- 
tions. The  initial  choice  of  E is  taken  to  be  E^;  thereafter,  the  new  is 
derived  from  the  old  as 

En+1  = En  + (-ED'EFn^ 

Once  the  elevation  E and  launch  speed  v^  are  determined,  the  orbit 
parameters  a,  e,  E^p,  AM  and  At  are  computed.  The  new  time  of  flight  is 
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compared  to  the  previous  value  and  the  whole  procedure  repeated,  unless  the 
difference  is  less  than  one  second. 

The  final  ECI  state  vector  is  then  transformed  to  the  launch  point 
"Radar"  XYZ  frame.  Launch  speed,  azimuth  and  elevation  are  computed  from  the 
Radar  XYZ  velocity  vector  to  provide  the  initialization  of  the  final  trajectory 
optimization  discussed  below  in  Section  6.3. 


6.2.3  Re-entry  Angle  Specified 

This  option  is  very  similar  to  th.  t described  above  for  the  launch  angle 
specified  case.  In  fact,  the  initialization  procedure  is  exactly  the  same, 
except  that  the  desired  launch  elevation  is  given  by 


ED  = I ErI  + °-017 


where  ER  is  the  re-entry  angle  in  radians  (ER  < 0).  The  additional  0.017 
radians  (1  degree)  is  the  result  of  observing  that  typically  re-entry  angles 
(magnitude)  are  smaller  than  the  launch  angles  by  0.01  to  0.02  radians. 

6.3  Launch  State  Vector  Optimization 


The  procedures  outlined  in  Section  6.2  above  lead  to  reasonably  good 
estimates  of  launch  state  vectors  based  on  a homogeneous,  spherical  earth 
model.  As  a consequence  of  the  simplified  earth  model,  the  target  trajectory 
will  not  impact  at  exactly  the  desired  coordinates.  The  problem  then  is  to 
perturb  the  launch  velocity  vector  in  such  a way  that  the  target  will  impact 
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at  the  desired  point  and,  in  addition,  will  satisfy  the  launch  angle,  re-entry 
angle  or  minimum  energy  constraint. 

An  iterative  technique  for  this  launch  state  vector  optimization  is 
presented  below  in  a general  mathematical  context.  Then  the  specific  cases  of 
interest- -launch  angle  specified,  re-entry  angle  specified,  minimum  energy-- 
are  considered  as  three  separate  special  examples. 


6.3.1  A Solution  of  Nonlinear  Equations 

The  impact  longitude,  impact  latitude  and  re-entry  angle  are  nonlinear 
functions  of  the  launch  velocity  vector  due  to  the  oblateness  of  the  earth  and 
due  to  the  gravitational  liaimonics  of  the  earth's  field.  In  this  section,  an 
iterative  technique  for  perturbing  the  launch  velocity  vector,  such  that  the 
trajectory  will  ultimately  impact  at  the  correct  latitude  and  longitude  with 
the  desired  re-entry  angle  (if  specified),  will  be  presented. 

Let  v be  the  launch  velocity  vector  and  let  y be  a vector  containing  the 
quantities- -impact  longitude,  impact  latitude  and  re-entry  angle.  Actually, 
it  is  the  vector,  q = y-y^,  which  is  to  be  minimized.  In  this  case,  y^  is  the 
vector  of  the  desired  values  of  y.  The  relation  between  q and  v can  be 
expressed  as 


q = f(v)  (6-44) 

where  f(v)  is  a nonlinear  function  of  the  vector  v.  The  objective  is  to  deter 
mine  v such  that  q = 0.  A perturbation  Av  of  the  vector  v results  in 
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q + A q s f(v+Av)  + JAv 


(6-45) 


where  J is  the  Jacobian 
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Now  let  vQ  be  an  initial  guess  (obtained  from  the  initialization  described  in 
Section  6.2).  Then, 


% ‘ f(v> 


(6-47) 


q - f(v  +Av  ) = ffv  ) + J Av 
H v o 0J  0/  0 0 


is  to  be  zero,  then  it  is  required  that 


f(0  + J_Av  = q + J Av  = 0 
0 0 0 o 0 


Av  = -J  q = -Ho 
o o no  o^o 


(6-48) 


-1 


where  H k J . Now  let 


h = VA 


(6-49) 


be  a new  estimate  of  the  Jacobian  where  A is  a perturbation  on  the  previous 
estimate.  Then 

<*i  = = VJiAuo 


or 


~ <?  +J«Av  +AAv 
“o  0 0 o 


h = AAV 


(6-50) 


since  J Av  = -J  H * -<?  . The  problem  now  is  to  determine  the  matrix  A 
o o o ono  no  r 

such  that  Eqn.  (6-50)  is  satisfied.  If  it  is  further  required  that 


AAu  = 0 


for  Au  orthogonal  to  Avq;  i.e., 


Av^Au  = 0, 


then  it  can  be  shown  that  A has  the  unique  solution 


T 

A0  .T 


(6-51) 


(6-52) 


Au* AV. 
o 0 


(6-53) 


This  gives  a new  value  for  the  Jacobian,  J 
required  is  the  inverse  of  Jx : 

1-1 


Jq+Aq,  but  what  is  actually 


Hi  - Ji 
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The  well  known  matrix  identity  (Householder  Identity) 


(I+uvV1  = I - uv 


l+vTu 


(6-55) 


permits  Hx  to  be  expressed  as 


Hi  = Ho 


Ho,  A v H 
on  1 oo 

X 


Av  Av  +Av„H  Oi 

0 o o o”  1 


(6-56) 


Thus,  when  initial  guesses  of  vq  and  Hq,  are  available  (or  computed) , the 
iterative  procedure  is  to  compute  q = f(vQ),  Avq  = -H0q0,  = vo+AvQ, 

etc.,  until  |q  | satisfies  the  desired  error  criterion.  The  initial  value  for 
the  vector  vQ  is  obtained  by  the  procedures  described  in  Section  6.2. 
Initialization  of  Hq  is  accomplished  by  perturbing  vQ,  one  component  at  a 
time,  and  computing  the  changes  in  the  components  of  q;  i.e., 
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Avok  = °»  W 


(6-57) 
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AV 


OJ 


The  resultant  matrix  J is  then  inverted  to  obtain  H . 

o o 

The  iterative  technique  for  the  solution  of  nonlinear  equations,  described 
above,  is  quite  general.  It  will  now  be  specialized  to  the  cases  of  interest. 

First  note  that  the  function  f(v)  used  to  relate  v and  q is  simply  the 
ballistic  equations  of  motion  described  in  Section  4.  The  velocity  vector,  v, 


Launch  speed 
Launch  azimuth 
Launch  elevation 


(6-58) 


is  converted  to  earth  surface  fixed  XYZ  coordinates  at  the  launch  point.  This 
is  followed  by  a conversion  to  an  ECI  state  vector  for  use  by  the  trajectory 
integration  routine,  lue  trajectory  is  next  integrated  4 either  the  re-entry 
altitude  or  the  impact  point  (or  both)  and  the  appropriate  values  for  the  y 
vector  are  computed,  where 


y = 


Re-entry  angle 
Impact  longitude 
Impact  latitude 


(6-59) 


The  q vector  is  simply  the  difference  between  the  computed  y vector  and  the 
desired  vector  y^. 
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6.3.2  Launch  Angle  Specified 


For  this  case,  the  launch  elevation  angle  is  held  fixed  at  the  desired 
value.  Launch  speed  and  launch  azimuth  in  an  earth  surface  (radar)  fixed 
frame  are  perturbed  as  described  above.  When  the  iterative  procedure  has 
converged  to  the  impact  point  (longitude  and  latitude) , the  launch  velocity 
vector  still  has  the  desired  elevation  angle  but  the  launch  speed  and  launch 
azimuth  are  altered  to  achieve  the  proper  impact  point  including  the  effects 
of  gravitational  harmonics  and  the  oblateness  of  the  earth  ellipsoid. 

The  constraint  on  re-entry  angle  does  not  apply  in  this  case,  so  the 
first  column  of  the  H matrix  is  set  to  zero  to  prevent  perturbations  of  the 
velocity  vector  v due  to  re-entry  angle  errors.  Also,  the  launch  elevation 
angle  is  not  included  in  the  perturbations  so  the  third  row  of  the  H matrix 
is  set  to  zero.  Thus, 


H - 0 


(6-60) 


It  is  readily  shown  that  t e method  of  updating  the  H matrix  (Eqn.  6-56) 
results  in  a propagation  of  the  null  elements  after  they  are  initially 
zeroed. 

Subroutine  NONLIN  (entry  point  SLV2)  is  called  to  perform  the  launch 
velocity  optimization  when  the  launch  elevation  angle  is  specified. 
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6.3.3  Re-entry  Angle  Specified 


For  this  case,  the  re-entry  angle  enters  the  computations  as  a constraint. 
Also,  the  launch  elevation  angle  is  perturbed  in  the  optimization.  The  H 
matrix  used  in  generating  launch  velocity  vector  perturbations  therefore  has 
all  non- zero  elements.  Convergence  of  the  iterative  procedure,  in  this  case, 
results  in  a launch  velocity  vector  which  impacts  at  the  desired  longitude 
and  latitude  and  which  has  the  requested  re-entry  angle. 

Subroutine  NQNLIN  (entry  point  SLV3)  is  called  to  perform  the  launch 
velocity  optimization  when  the  re-entry  angle  is  specified. 

6.3.4  Minimum  Energy  Optimization 

The  procedure  for  launch  velocity  vector  optimization  subject  to  the 
minimum  energy  constraint  is  somewhat  more  complicated  than  the  above  two 
cases.  It  is  first  assumed  that  the  initialization  described  in  Section  6.2.1 
is  close  to  the  true  minimum.  Then,  the  launch  angle  specified  option  is 
called  three  times  with  E = E . , E = E ,,+AE  and  E = E^+-AE,  where  E . is 
the  launch  elevation  angle  derived  from  the  initialization  described  in 
Section  6.2.1  and  AE  is  currently  set  to  0.005  radians.  The  three  launch 
speed-launch  elevation  pairs  are  next  fitted  with  a parabola  to  determine  the 
elevation  which  minimizes  the  launch  speed.  A final  call  to  Subroutine  NONLIN 
(entry  point  SLV2)  with  this  elevation  angle  results  in  a launch  velocity 
vector  corresponding  to  a minimum  energy  trajectory  between  the  specified 
launch  and  impact  points. 
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6.4  Trajectory  Information 


The  presentation  given  above  concerns  the  generation  of  target  trajectories 
subject  to  certain  constraints.  There  is  the  alternate  possibility  that  a 
target  state  vector  for  a particular  trajectory  is  available  for  input  to  the 
program  directly.  In  this  case,  it  is  desirable  to  determine  various  parameters 
(e.g. , launch  point,  impact  point,  re-entry  angle)  of  the  trajectory  for 
documentation  of  a particular  run  of  the  ARIES  Program.  The  mathematics  for 
computation  of  the  relevant  trajectory  parameters  are  presented  below. 

6.4.1  Input  State  Vectors 

There  are  eight  options  for  specification  of  a target  trajectory  as 
explained  in  Reference  1,  Section  S.3.  Options  6,  7 and  8 are  the  minimum 
energy,  launch  angle  specified,  and  re-entry  angle  specified  options  which 
require  the  trajectory  generation  processing  discussed  above.  The  other  five 
options  require  the  input  of  a complete  state  vector  in  a particular  coordinate 
frame  as  follows. 

1.  EC I XYZ 

2 . Radar  XYZ 

3 . Radar  RAE 

4.  Phased  array  XRF 

5.  Phased  array  RUV 

Each  of  these  five  options  requires  a state  vector  time  of  validity  (tag  time) . 
Options  2 through  5 also  require  the  radar  longitude,  geodetic  latitude  and 
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height  above  the  earth  ellipsoid  to  be  specified  in  the  input  stream.  Options 
4 and  5,  in  addition,  require  the  azimuth  and  elevation  angles  of  the  phased 
array  boresight.  The  ECI  State  Vector  (Option  1)  is  referenced  to  a particular 
longitude  at  the  launch  time  (TAL  ■ 0);  this  longitude  is  a required  input. 
(Option  9 for  a satellite  orbit  has  not  been  implemented.) 

The  first  step  in  the  handling  of  these  input  state  vectors  is  to  call  a 
series  of  coordinate  transformation  subroutines  to  finally  obtain  a state 
vector  in  ECI  coordinates  referenced  to  the  Greenwich  meridian.  This  ECI 
State  Vector  is  then  integrated  backward  in  time  to  the  launch  point  (Sub- 
routine LAUNCH)  and  forward  in  time  to  the  altitude  of  re-entry  and  to  the 
impact  point  (Subroutine  IMPACT) . The  State  Vector  returned  from  Subroutine 
LAUNCH  has  a tag  time  referenced  to  time  of  Hunch  (TAL). 


6.4.2  Launch,  Impact  and  Re-entry  Parameters 


Both  the  LAUNCH  and  the  IMPACT  Subroutines  require  trajectory  integration 
to  a specified  altitude  (zero  for  launch  and  impact  points,  300  Kft.  for 
re-entry).  At  the  specified  altitude,  the  state  vector  is  appropriately 
transformed  to  determine  the  target  longitude  and  geodetic  latitude.  The 
launch  elevation  angle,  impact  angle  or  re-entry  angle  is  also  determined  as 
the  angle  between  the  velocity  vector  and  the  local  horizontal  plane.  The 
azimuth  angle  of  the  velocity  vector  at  launch  is  also  determined.  Each  of 
these  computations  is  performed  in  a radar  XYZ  frame  centered  at  the  longitude, 
latitude  and  altitude  of  the  target. 


Consider  an  ECI  target  state  vector,  Xj,  somewhere  along  a trajectory. 
The  altitude  of  the  target  (above  the  earth  ellipsoid)  is  given  approximately 
by  (Eqn.  2-14) : 


h = r 


(0-61 


where  r =yx|+y|+z|,  cosA  = p/r  and  p = ^Xj+y|\  (As  before  is  the  North 
polar  radius  and  e is  the  eccentricity  of  the  earth  ellipsoid.)  The  angle  A 
is  approximately  equal  to  the  geocentric  latitude  of  the  point  on  the  surface 
of  the  earth  from  which  altitude  is  measured  (i.e.,  the  normal  to  the  earth 
ellipsoid  passing  through  the  target  position) . A first  approximation  to 
the  geodetic  latitude  is  therefore 


<t>  - tan 


The  longitude  of  the  target  is  given  by 


r zi  i 

(6-62)  ; J 
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X = tan  — - u it,  Os  is  2tt 
/ e 


(6-63) 


where  the  second  term  accounts  for  the  rotation  of  the  earth  from  the  reference 
time  of  the  ECI  frame  (usually  TAL  = 0)  to  the  tag  time  of  the  state  vector. 

The  parameters  X,  $ and  h give  the  position  of  the  target  with  respect  to 
the  surface  of  the  earth.  It  is  also  desirable  to  determine  the  azimuth  and 
elevation  of  the  velocity  vector  with  respect  to  an  earth  surface  fixed  frame 
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at  this  position.  To  acconplish  this,  the  following  three  unit  vectors  are 
defined: 


lI  - 


*1  , 


*1  - cos<v  + r M + sin*  s 


(6-64) 


yi  . xi  . 

*2-  -r^x  + r^y 


(6-65) 


. /XI  . yI  • \ 

*3  = -s^(-  *x  + - + cos* 


(6-66) 


The  unit  vector  x.,  is  normal  to  the  earth  ellipsoid  at  latitude  <t> ; i2  is  normal 
to  -t1 , parallel  to  the  equatorial  plane  and  points  to  the  east ; i3  is  normal  to 
both  i-x  and  -l2  and  points  to  the  north.  That  is;  ix , i2  and  t3  correspond  to 
the  Z,  X and  Y axes  of  a radar  XYZ  frame.  The  conponents  of  the  velocity 
vector  in  the  ECF  frame  are 


% = Vvi 

S’f  = >V"exI 


iF-  il 


(6-67) 


The  target  speed  in  the  ECF  frame  (or  radar  XYZ  frame)  is 


v =/xF+yF+*F’ 


(6-68) 
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The  components  of  the  velocity  vector  in  the  earth  surface  fixed  frame  are 
obtained  by  taking  the  following  dot  products 


vz  = = cos<() 


piWiI 


+ sin<{>  t. 


(6-69) 


Vx  = VyF  = 


■wVf 


(6-70) 


vy  = 


. [V 

P = -sinct>  I — 


I+W 


+ COS<()  t. 


(6-71) 


where  Vp  = •tx*p+<C^p+^z2p‘  Note  that  the  effects  of  the  rotation  of  the  earth 

cancel  in  the  expressions  for  v and  v . The  elevation  angle  of  the  velocity 

z y 

vector  (which  may  be  the  launch  angle,  impact  angle  or  re-entry  angle)  is  given 


p • -ir* 

Ev  = sin  - 


(6-72) 


and  the  azimuth  angle  of  the  velocity  vector  is  given  by 


-i/vx 


\ = tan 


(6-73) 


The  computation  of  A,  <)>,  h,  E and  Ay  for  an  ECI  state  vector  are  per- 
formed in  Subroutine  ANGLES.  Subroutines  LAUNCH  and  IMPACT  both  call  Subrou- 
tine ANGLES  after  the  ECI  state  vector  has  been  integrated  forward  or  backward 


to  the  relevant  target  altitude.  Various  parameters  of  launch  (\,<|>,E  ,Ay) 
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re-entry  (X,4>,EV)  and  impact  (x,<j>,Ev)  are  saved  in  the  target  information 
files. 

Section  2 of  the  ARIES  Test  Report  for  a particular  run  of  the  Program 
(see  Reference  1,  Appendix  II)  contains  trajectory  information  for  each  target, 
In  particular,  the  following  information  is  available  in  the  Test  Report: 


1.  Initial  state  vector  data 

2.  Initial  conditions  for  target  tumbling  and  RCS  computations 

3.  Launch  parameters 

4.  Re-entry  parameters 

5 . Impact  parameters 


SECTION  7- -MAXIMUM  LIKELIHOOD  TRAJECTORY  ESTIMATION 

At  present  the  only  "tracking"  algorithm  built  into  ARIES  is  the  Maximum- 
Likelihood  Estimator  (MLE) . The  MLE  does  not  provide  estimates  of  a target's 
position  after  each  radar  measurement,  as  do  recursive-type  tracking  schemes 
such  as  the  a -6  tracker  or  a Kalman  filter.  Instead,  the  MLE  provides  ARIES 
only  with  an  endpoint  state  vector  estimate,  valid  at  the  end  of  the  track 
interval (s)  of  interest.  During  the  "track"  interval  itself,  the  MLE  Subrou- 
tine MAXLIK  simply  stores  up  the  (noisy,  biased)  radar  position  measurements, 
along  with  the  variances  01  each  measurement. 

Once  all  such  measurements  have  been  accumulated,  the  MLE  attempts  to 
find  that  endpoint  state  vector  estimate  which  minimizes  the  "residual  errors" 
between  the  physical  trajectory  (uniquely  defined  by  the  endpoint  state  vector) 
and  the  set  of  radar  measurements.  The  error  used  in  this  minimization  is 
formed  as  the  sum  of  the  weighted,  squared  differences  between  the  physical 
trajectory  and  the  set  of  radar  measurements.  Each  error  residual  is  weighted 
according  to  the  standard  deviation  of  the  corresponding  position  measurement. 
Since  each  error  residual  is  weighted  according  to  the  variance  of  the 
| corresponding  position  measurement,  the  final  state  vector  estimate  corresponds 

| . to  a minimum  variance  estimate.  This  minimum  variance  estimate  is  also  a 

: maximum  likelihood  estimate  of  the  trajectory,  when  the  measurement  errors 

have  a zero-mean,  Gaussian  distribution.  If  there  are  no  bias  errors  in  the 

data  (such  as  those  due  to  inadequate  calibration  or  those  due  to  imperfect 

I 

i correction  of  tropospheric  or  ionospheric  refraction  effects),  the  trajectory 
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estimation  procedure  described  below  corresponds  to  Maximum  Likelihood 
Estimation.  An  extension  of  the  procedure  to  include  estimation  of  biases 
would  again  result  in  Maximum  Likelihood  Estimation  of  both  the  state  vector 
and  the  biases. 

7.1  The  MLE  Method 

Let  xg  be  a reference  ECI  state  vector  at  an  arbitrary  reference  time  tg. 
(Note  that  xg  is  the  state  vector  to  be  perturbed  to  achieve  the  "optimal", 
weighted  least- squares  trajectory  fit  to  the  measurement  data.  For  this  study 
the  state  vector  is  valid  at  the  endpoint  of  the  tracking  interval,  but  this 
is  not  a basic  constraint.)  Also  let  the  set  of  radar  measurements  at  N 
points  in  time  t , t2,...tN,  be  denoted  by  the  vectors  -^Ct^),  i=l,  2,...N. 
These  radar  position  measurements  may  be  either  in  RAE  or  RUV  coordinates 
dependent  on  whether  the  radar  is  a dish  type  or  a phased  array  type.  The 
reference  ECI  state  vector  xs  is  integrated  successively  to  the  times  of  the 
radar  measurements.  At  each  of  these  times,  the  state  vector  xs(t^)  is 
transformed  to  the  appropriate  radar  frame  (RAE  or  RUV)  to  form  the  set  of 
"estimated"  position  vectors,  *s(tp.  It  is  the  weighted,  squared  differences 
between  \(t^)  and  ^s(tp  which  are  minimized  in  the  trajectory  fitting 
process. 

It  has  been  shown  in  Section  5.2  that  the  transformation  from  an  ECI  state 
vector  into  radar  XYZ  (East-North-Up)  coord iuates  can  be  performed  via: 


xsr(t)  = A(<M)G(t)xs(t)+b(<fr) 


(7-1) 
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where  the  transformation  matrix  G(t)  accounts  for  the  rotation  of  the  earth  in 
time  t (see  Eqn.  5-1).  The  transformation  matrix  A(4>,X)  and  the  translation 
vector  6(<j>)  are  those  used  to  transform  from  an  ECF  frame  to  a radar  XYZ  frame 
(see  Eqns.  5-14  and  5.15).  For  a phased  array  radar,  the  additional  trans- 
formation from  radar  XYZ  to  phased  array  XYZ  (see  Eqn.  5-26)  must  be  included 
so  that  the  vectors  xgr  are  defined  in  the  proper  coordinate  system  for  the 
radars . 

Also,  we  have  denoted 

(7-2) 


where  the  time  dependence  is  implicit  in  these  definitions. 

From  Equation  (7-1)  it  is  clear  that  a perturbation  in  xg  will  lead  to  a 
corresponding  perturbation  in  xgr: 

5xsr  = AG  6xs  (7-4) 

Actually  we  want  to  relate  perturbations  in  the  reference  state  vector  xg  to 
perturbations  in  the  "estimated"  target  position  n . To  do  this  for  the  RAE 
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case,  we  use  Eqn.  (5-21)  in  Eqn.  (5-20)  to  determine  the  relation  between 
perturbations  of  *g  and  xgr: 


6*s  = C 6xsr 


(7-5) 


where  the  matrix  C is  given  by 


C ft 


cosEsinA 

cosA/RcosE 

-sinEsinA/R 


cosEcosA  sinE 
-sinA/RcosE  0 
-sinF,cosA/R  cosE/R 


(7-6) 


For  a phased  array  radar  (RUV  coordinates)  the  matrix  C is 


RU 

RV 

RW 

1 

C ft  ~ 
L R 

1-U2 

-UV 

-uw 

(7-7)  | 

-uv 

1-V2 

-vw 

where  W = ^Tu2-^1.  Note  tha^  the  matrix  C varies  as  the  target  moves  along 
its  trajectory;  for  simplicity,  this  time  dependence  has  not  been  explicitly 
included  in  the  above  relations. 

Combining  Equations  (7-4)  and  (7-5)  one  can  then  relate  perturbations  in 
the  position-measurement  vector  directly  to  perturbations  in  the  ECI  state 
vector: 


6*5  = CAG«xs  (7-8) 

Back  to  the  NILE  problem:  Given  the  reference  state  vector  xg,  the  MLE  Sub- 

routine MAXLIK  first  uses  the  equations  of  motion  (Section  4.3)  to  obtain  xg 





72 


and  then  Equation  (7-1)  to  obtain  xgr  at  each  measurement  time.  Next  the  *sr 
state  vector  is  transformed  into  a (the  set  of  target  position  estimates 
based  upon  the  ECI  reference  vector  xs) , using  the  appropriate  transformations 
(given  in  Section  5.3  for  an  RAE  radar,  and  in  Sections  5.4  and  5.5  for  an  RUV 
radar) . 

In  deriving  the  MLE  estimator,  we  assume  that  perturbations  in  the  refer- 
ence state  vector  positions  x and  rates  x at  time  t cause  corresponding 

o i)  5 

perturbations  in  xs(t.)  via  the  equation 

«xs(ti)  = Sxs  + (ti-ts)«*s  (7'9) 

This  assumption  becomes  increasingly  valid  as  6xg  ard  6*s  become  small,  as 
they  do  on  successive  MLE  iterations. 

The  result  of  the  above  discussion  is  as  follows:  Given  perturbations 

about  the  end-of-track  ECI  position  estimate  xr  and  rate  estimate  i at  a time 
ts,  we  may  determine  the  corresponding  perturbation  in  the  target  position 
estimate  at  time  t^  in  radar  measurement  coordinates  using  Equations  (7-8)  and 
(7-9)  to  obtain: 

<5/ts(ti,6xs,6xs)  = F(ti)[6xg  + (ti-ts)6^s]  (7-10) 

where  we  have  defined 

F(tp  = C(ti)A  G(t.)  (7-11) 

At  this  point  we  have  a set  of  radar  measurement  data  /i  (tp  , a set  of  "esti- 
mated" target  positions  K (t-)  based  on  the  reference  ECI  state  vector  x.(t.) 

3 1 s s 
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is: 


I: 


!■ 


integrated  to  the  appropriate  times,  and  a set  of  position  error  estimates 


^sCt;)  from  Eqn.  (7-10).  These  are  now  combined  in  a quadratic  form  as 


follows: 


N 


Q - L^VV  " W • 6/ts(ti»6xs»6y } 


- «*s(ti,«Xg,«*s)) 


(7-12) 


where  W is  a weighting  matrix  taking  into  account  the  quality  of  the  radar's 
position  i,.easurements.  Thus,  for  an  RAE  radar: 


OKt.)  = 


0 


0 


°R(ti) 


^(tf) 


«i<V 


(7-13) 


where  a^(t^)  is  the  range  measurement  variance,  aj^(t,. ) is  the  azimuth  measure- 
ment variance.,  and  a£(t^)  is  the  elevation  measurement  variance  for  the  radar 
measurements  at  tune  t^.  For  a phased  array  radar  (RUV  data),  is  replaced 
by  a2  and  o£  by  a2.  For  simplicity,  let  the  arguments  t-  be  simply  denoted 
by  the  index  i,  Then  with  the  additional  definitions 
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"M 


M 


I 


m 


i 


■m 


•m 


s 


n 

A 


«*(i)  = %(i)  - Ag(i), 


(7-14) 


i^(i,(Sxs,6xs)  = 5A(i)-6/ls(i,6xs>6xs), 


(7-15) 


Eqn.  (7-12)  becomes: 


Q = £ j((i,<5xs,6xs)  IW(i)f^(i><Sxs,6xs)] 


(7-16) 


If  the  components  of  6xg  are  denoted  by  6xg>  6yg,  6zs,  then  to  minimize  Q one 
must  set  the  partial  derivatives  of  Q with  respect  to  6x_,  6y  and  6z  to 


zero: 


an  T ^ (i  > 

I?  - 2El((i,*<s.*y™  — - 


(7-17) 


here  a takes  on  the  values  6xs,  6yg  and  <5zs>  Three  simultaneous  equations 
are  generated  from  Eqn.  (7-17)  which  may  be  combined  into  the  single  matrix 
equation 


£U(i,o*  s,6xJ]TW(i)F(i)  = (0,0,0) 
i=l  * s 


(7-18) 
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Similarly,  if  the  components  of  6x  , are  denoted  by  , 6y  , Si  , and  if  the 

to  6fs  and 

the  following  matrix  equation  results: 


partial  derivatives  of  Q with  respect  to  6*s,  and  62s  are  set  to  zero,  then 


N 


£[(5(i,6xs,6is)]Ta;(i)F(i)(ti-ts)  = co,o,o) 


(7-19) 


The  ML  Estimator  simply  has  to  solve  these  linear  equations  for  (xs  and  6*s; 
then  the  end-of-track  state  vector  estimate  can  be  updated: 


xs,new  xs,old  + "xs 


(7-20) 


xs,new  *s,old  + 6xs 


(7-21) 


The  matrix  equations  for  the  ML  Estimator,  Eqns.  (7-18)  and  (7-19),  can  be 
written  in  a more  compact  form  as  follows.  First  take  the  transpose  of  these 
equations  and  substitute  the  expression  for  |{(i,6xs,6*s)  from  Eqn.  (7-15). 

Then  define  the  following  vectors  and  matrices 


rf,  - E 


i=l 


FT(iMi)6*(i) 


(7-22) 


N 


d?  = E(trts)FT(i)W(i)6^(i) 
i=l  1 5 


(7-23) 


These  are  not  to  be  confused  with  the  A and  b which  were  defined  above  in 
Equation  (7-1). 

Note  that  the  ML  Estimator  is  an  iterative  procedure:  An  initial, 
reference  state  vector  xg  is  determined  (equal  to  the  true  target  state  vector 
in  ARIES) , then  perturbations  Sxg  and  6xg  are  determined  by  the  above  pro- 
cedure. The  new  value  of  from  Eqns.  (7-20)  and  (7-21)  is  then  used  to 
repeat  the  process.  Since  ARIES  starts  with  an  excellent  value  for  xs>  the 
convergence  of  this  process  is  very  rapid;  in  fact,  only  two  iterations  are 
used. 

7.2  Error  Residuals 

The  residual  errors  between  the  target-position  measurement  for  the  radar 
*.m(i)  and  the  target-position  prediction  *s(i)  at  each  time  ti  (where  *s(i)  is 
computed  from  xg,  the  final -iteration  value  of  the  MLE  end-of-track  state 
vector  estimate)  are  computed  via 

NFP1AX 

Mean  square  range  residuals  = £ [ 6r . (i) ] 2 (7-32) 

i=l 

where  NPTMAX  is  the  number  of  measurements  made  during  the  Monte  Carlo  run, 
and  where  6rx  is  the  first  component  of  the  residual  error  vector  6A.  defined 
in  Equation  (7-14).  Similar  mean  square  residuals  are  computed  for  the  other 
two  components  of  £*.  These  mean  square  residuals  are  computed  for  each  Monte 
Carlo  run,  and  then  averaged  over  all  runs.  The  resultant  mean- square -resid- 
uals are  printed  out  in  Section  11  of  the  ARIES  Test  Report. 


78 


SECTION  8- -HANDOVER  STATE  VECTOR  QUALITY 


In  ARIES,  the  MLE  state  vector  estimate  of  the  position  and  velocity  of  a 

target  at  the  end  of  a tracking  interval  (or  extrapolated  ahead  in  time  beyond 

the  end  of  the  tracking  interval)  will  not,  of  course,  correspond  exactly  to 

the  true  position  and  velocity  of  the  target,  due  to  the  various  noises  and 

biases  which  were  present  in  the  radar  measurements,  as  explained  fully  in 

Reference  2.  The  error  vector  at  any  time  between  the  extrapolated  state- 

vector  estimate  x and  the  true  state  vector  x.  is  defined  as 

true 


e - v - x. 


(8-1) 


where  both  x and  *true  are  6 -component  ECI  state  vectors. 

The  volume  spanned  by  the  error  vector  £ is  of  interest  to  the  command 
guided  intercept  problem,  since  it  represents  the  uncertainty  in  the  handover 
intercept  position  of  the  target.  If  the  interceptor  is  being  guided  by  a 
different  radar  than  the  one  which  tracked  the  target  (i.e.,  a radar  with  a 
different  set  of  biases,  tracking  at  a different  elevation  with  different 
amounts  of  refraction,  etc.),  then  the  error  as  defined  in  Equation  (8-1)  is 
of  chief  interest,  as  described  below  in  Section  8.1.  If,  however,  the  same 
radar  is  tracking  both  the  target  and  the  interceptor,  then  the  biases  would 
tend  to  wash  out,  and  the  uncertainty  after  subtracting  out  the  mean  error  is 
of  more  significance,  as  described  in  Section  8.2. 
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8.1  Correlation  Matrices 

On  the  mth  Monte  Carlo  simulation  run  the  ARIES  Test  Program  computes  a 
"correlation"  matrix  estimate  via 


m 


(8-2) 


where  is  the  error  defined  above  in  Equation  (8-1)  on  the  mth  Monte  Carlo 
run.  These  correlation  matrix  estimates  are  then  averaged  over  all  runs: 


c A 


MC 

NT  £ 
m=l 


1 


m 


(8-3) 


It  is  clear  from  Equation  (8-2)  that  each  of  the  correlation  matrices  Cm 

is  symmetric;  hence  C is  symmetric  and  only  the  "lower  half"  is  printed  out 

in  Section  12  of  the  ARIES  Test  Report,  as  explained  in  Reference  1. 

For  the  error  ellipsoid  calculations  (see  Section  3.3),  it  is  necessary 

to  diagonalize  the  "position"  correlation  matrix  P by  means  of  an  orthogonal 

transformation  (P  is  the  upper  left-hand  quadrant  of  C;  P is  a 3 x 3 matrix). 

This  diagonal izat ion  is  only  possible  if  the  matrix  P is  positive  definite; 

T 

i.e.,  if  the  quadratic  form  x Px  is  non-negative  for  all  real  values  of  the 
variables  x^,  and  is  zero  only  if  each  of  the  x^  variables  is  zero.  For  the 
P matrix,  this  statement  translates  to  the  condition 


(x 


. e,  +x„e  +x„e.,„)2 
i im  2 2m  3 3m' 
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which  is  obviously  non-negative.  If  MC=1  or  MC=2,  it  can  be  readily  shown 
that  there  are  an  infinite  number  of  non-zero  values  of  the  x^'s  for  which  the 
above  quadratic  form  is  zero.  This  violates  the  conditions  for  P to  be  posi- 
tive definite.  For  MC  > 3,  it  is  not  generally  possible  to  make  the  quadratic 
form  equal  to  zero  with  non-zero  values  for  the  x^'s.  Exceptions  are  possible 
if  the  error  vectors  are  linearly  related;  this  is  unlikely  in  the  case  at 
hand,  particularly  for  large  values  of  MC.  The  reason  why  the  position 
correlation  matrix  is  singular  for  one  or  two  Monte  Carlo  runs  can  be  simply 
explained  as  follows:  a volume  in  space  requires  three  non-parallel  vectors 

for  its  definition  (one  error  vector  defines  a point,  two  error  vectors  define 
a plane) . Consequently,  for  the  ARIES  Program  to  generate  meaningful  error 
volumes  based  on  position  errors,  it  is  necessary  to  perform  a minimum  of 
three  Monte  Carlo  runs. 


8 . 2 Covariance  Matrices 

If  the  mean  error  z (obtained  by  averaging  over  all  Monte  Carlo  runs)  is 
subtracted  from  the  error,  an  estimated  covariance  matrix  could  be  obtained  on 
each  Mon+  3 Carlo  run  via 


Vm  4 <Ve’ 


(8-4) 


and  then  averaged  over  all  Monte  Carlo  runs  to  obtain 


- A , MC  A 

m='. 


(8-5) 
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In  practice,  what  is  actually  done  by  ARIES  is  to  first  compute  C via  Equation 
(8-3),  and  next  compute 


1 MC 

e = MC  S,  % 

m=l 


(8-6) 


and  finally  to  obtain  V using  the  relation 


V = C - e eT 


(8-7) 


As  with  the  correlation  matrix  C this  covariance  matrix  V is  symmetric  and 
thus  only  the  lower  half  appears  in  Section  12  of  the  ARIES  Test  Report, 

The  "position"  covariance  matrix  P (the  upper  left-hand  3x3  sub-matrix 

A 

of  V)  is  singular  for  MC  < 4;  the  process  of  subtracting  the  means  of  the 
errors  effectively  reduces  the  dimensionality  of  the  space  by  one.  Conse- 
quently, for  the  ARIES  program  to  generate  meaningful  error  volumes  based  on 
the  position  error  covariance  matrix,  it  is  necessary  to  perform  a minimum  of 
four  Monte  Carlo  runs. 


8.3  Handover  Error  Ellipsoid 

The  correlation  and  covariance  matrices  described  in  the  preceding 
sections  are  for  the  errors  as  measured  in  the  basic  ARIES  ECI  coordinate 
system  described  earlier  in  Section  3.1.  The  drawback  of  such  matrices  is 
that  all  elements  are  non- zero  - that  is  the  X,  Y,  7,  error  components  are 
correlated  in  this  ECI  frame. 


,.J  - . .--.it!  >•  ■-■•Ai-A'ite t/A*  I.*,  v,  , „ , 


To  eliminate  the  correlations  between  the  error  vector  components  requires 
the  definition  of  a new  coordinate  frame.  Let  a "position”  correlation  matrix 


P be  defined  as  the  upper  left-hand  quadrant  of  a correlation  matrix  C (or  of 


a covariance  matrix  V ) 


[r] 


(8-8) 


Then  we  want  to  determine  a new  coordinate  system  such  that  the  ECI  position 

T T 

error  vector  e = (e^e^e^)  is  related  to  the  position  error  vector  (e')  = 

(ej.e^eg)  in  the  new  coordinate  system  via  the  transformation  matrix  Q.: 


e'  = Qjz 


(8-9) 


-1  T 

where  Q is  the  normalized  modal  matrix  of  P with  the  property  that  Q = o . 
The  new  "position"  correlation  matrix  P*  is  then  easily  shown  to  be  related  to 


P via 


P'  = dTPd 


(8-10) 


Given  the  symmetric,  nonsingular  matrix  P,  the  problem  is  to  find  the 
matrix  Q.  such  that  Q.  PQ  is  a diagonal  matrix.  This  is  a standard  problem 
(c  f,,  Hildebrand,  Reference  5),  equivalent  t.o  determining  the  eigenvalues  of 
the  matrix  P via  setting  the  following  determinant  to  zero: 


P-U  ■ 0 


(8-11) 
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where  A is  an  eigenvalue.  For  this  case  there  are  three  such  eigenvalues; 
they  are  the  roots  of  the  cubic  equation 


1 l'x 

P1 2 

P13 

21 

P22‘X 

P2  3 

= -A3-pA2-qA-r  = 0 

(8-12) 

31 

P32 

P -A 

33 

where 


p*  -(pn+p22*p33> 

^ = P22P33+P33P11+P11P22'P23‘P13'P12 


and 


p2  p +p2  p +p2  p _2P  P P -P  P P 
13  22  12  33  23  11  12  23  13  11  22  33 


The  solution  of  (8-12)  is  well  known.  Define 


(8-13) 

(8-14) 


(8-15) 


a = q - p2/3 

b = r-§a.‘7p3 

For  a positive  definite  matrix,  compute 


and  obtain  the  eigenvalues  as 


(8-16) 

(8-17) 


(8-18) 


(8-19) 
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s? --*  ts-- ^^jajfrWT^ .'.vy,;  ^v --y:. at^p. ,-. c *:?j27 -.r.v.-A-y  : r * - -.-y : 


cos  1 7 + r)  I - ? 


-[2f? 

= [2/?cos (j-  r)  "I 


(8-20) 

(8-21) 


Since  the  matrix  P is  positive  definite,  the  eigenvalues  (X^X,^,)  are  all 
real  and  positive.  It  then  follows  that  P1  has  the  desired  diagonal  property 


P = 


\l  0 
0 X 
0 


0 

0 


2 

0 X 


3J 


(8-22) 


The  square  roots  of  the  diagonal  elements  of  such  diagonalized  "position" 
correlation  (or  covariance)  matrices  can  be  thought  of  as  the  semi-axes  of  an 
ellipsoidal  error  volume  - that  is,  the  errors  in  this  new  coordinate  frame 
are  uncorrelated  and  have  rms  values  /7  and  X} , X0  and  X3  are  also 

printed  out  in  Section  12  of  the  ARIES  Test  Report. 

8.4  Handover  Error  Sphere 

While  the  coordinate  transformation  Q used  in  the  preceding  section 
rendered  uncorrelated  errors,  the  eigenvalues  x X2  and  X3  alone  are  insuffi- 
cient to  specify  the  uncertainty  volume  of  the  target  state  vector  at  hand 
over.  The  orientation  of  the  ellipsoid  axes  must  also  be  taken  into  account. 
In  order  to  avoid  this  complexity,  it  is  often  more  useful  to  think  of  an 
"uncertainty  sphere"  - namely  a sphere  of  radium  R,  centered  at  the 
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(extrapolated)  handover  point,  with  R chosen  such  that  the  error  vector  lies 
within  the  sphere  with  probability  = 0.99. 

The  diagonalization  of  the  position  correlation  (or  covariance)  matrix 
assures  us  that  the  error  components  in  the  new  coordinate  frame  (see  preceding 
section)  are  uncorrelated.  If  further  we  assume  these  errors  to  have  "normal" 
or  "Gaussian"  distributions,  then  the  probability  of  the  error  vector  e ' 
(defined  in  Equation  (8-9)]  is  given  by 


P(e')  s 


X2 

IxT 


Y2 

^2 


z2 

2'x  3 


(2tt)1*5v^7 


^2  V? 


A P(X,Y,Z) 


(8-23) 


(where  for  notational  convenience  the  components  ej,  e'2 , e^  of  the  error 
vector  e'  have  been  replaced  by  X,  Y and  Z,  respectively) . The  probability  of 
the  error  vector  falling  within  a sphere  of  radius  R is 


R /r2-X2'  V^2-X2-Y2’ 


PCI e*  | < R) 


' 8/dx/  dY  / 


p(X,Y,Z)dZ 


(8-24) 


where  due  to  the  symmetry  we  only  have  to  integrate  over  one-eighth  of 
the  sphere.  If  we  use  Equation  (8-23)  in  (8-24),  along  with  a set  of 
normalized  variables  defined  by 


Jn? 


( 8-26) 


\FT 


(8-27) 


then  we  obtain 


P(k’|  < R) 


L 


h-^rf-V-$C> 

dy  e~Y  I e'7-2  dz 


(8-28) 


where  we  have  defined: 


cx  - R/Jn? 


Cy  - RA/^ 


(8-29) 


(8-30) 


Cz  " R/>^7 


(8-31) 


Note  that  the  integration  is  now  over  an  ellipsoid  rather  than  a sphere, 
If  we  transform  y and  z to  polar  coordinates  r and  0,  where 


-v'F? 


(8-32) 


0 = tan  (z/y) 


(8-33) 


then  we  may  replace  the  two  inner  integrals  of  Equation  (8-28)  with 
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(8-35) 


Use  of  the  double  angle  formula,  2sin20  = l-cos20,  in  Eqn.  (8-35)  results  in 
the  following  relation 


(8-36) 


or, 


where 


r ^ F1+k2cos20 


(8-37) 


(8-38) 
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VrjTi/i'T.^&V;;^'.^*:*^  ;;  - - -.  ***•*'  ' 


and 


k = 
k2 


ry 


2,1  - kry 


(8-39) 


Substitution  of  Eqn.  (8-37)  into  Eqn.  (8-34)  gives 


ti/2  r(0) 


/*/ 


-r2,  1 

re  dr  = 4 


TT 


17  J | k1+k2cos(|) 


d<J> 


(8-40) 


*-  0 


where  the  change  of  variables  20  = <t>  has  been  used. 

When  Equation  (8-40)  is  used  to  replace  the  inner  two  integrals  of  (8-28), 
we  obtain 


PCI 


< R) 


- TTs  / *<  *‘X"  | * -Jex pj 


-1 


k1+k2cos"i}) 


d<(>' 


(8-41) 


The  double  integral  in  Eqn.  (8-41)  cannot  be  reduced  any  further;  there- 
fore, it  is  evaluated  numerically  by  a Gauss  quadrature  method. 

If  we  specify  an  error  sphere  radius  R,  Equation  (8-41)  allows  the  compu- 
tation of  the  probability  that  the  error  vector  lies  within  the  sphere.  How- 
ever, it  is  usually  desired  to  find  the  value  of  R which  gives  a probability 
of  0.99.  This  value  of  R is  found  numerically  by  means  of  the  secant  iterative 
method  described  in  Reference  4. 
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